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Summary 


..iiis  report  summarizes  the  recent  work  in  the  application  of 
the  LIE-series  method  to  the  solution  of  ordinary  and  partial 
differential  equations. 

After  a  short  introduction  the  pov;er  scries  method  which  is  a 
special  case  of  the  Lie  series  method  of  chapter  III  is  de¬ 
scribed  in  chapter  II.  Further  v;c  discuss  the  interesting  con¬ 
cept  cf  recursion  formulas  and  the  calculation  of  the  "trans¬ 
fer  matrix"  (connection  matrix),  the  derivatives  of  the  solu¬ 
tion  with  respect  to  the  initial  values. 

Chapter  III  deals  with  the  numerical  evaluation  of  the  Lie 
series  perturbation  formula.  This  chapter  contains  the  result r 
of  the  report  /29/,  which  has  been  v/ritten  together  with  II. 
Knapp  at  the  MRC,  Madison,  Wisconsin.  Suitable  quadrature  for¬ 
mulas  and  recursions,  statements  on  the  order  and  error  esti¬ 
mation  are  given.  Numerical  examples  finish  the  chapter  and 
compare  the  method  also  with  that  of  Fehlberg. 

•n  chapter  IV  wo  prove  Ortibner’s  integral  equation  which  leadr 
to  short  proofs  of  the  formulas  of  chapter  III  and  to  various 
generalizations  of  the  method.  A  survey  of  those  is  presented 
at  the  end  of  this  summary. 

Chapter  V  Generalizes  the  concept  of  Runge-Kutta  to  methods 
with  multiple  nodes,  which  is  possible  with  the  use  cf  the  Lie 
differential  operator  D.  A  general  theory  is  developed  and  the 
method  of  Fohlberg  is  shown  to  be  a  special  case. 

Chapter  VI  dealQ  with  the  stop-size  control  and  chanter  VII 
shows  the  application  of  centralized  Lie  scries  to  the  calcu¬ 
lation  of  switch-on  transients  occuring  in  the  telegraphic 
equation. 
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Chapter  I  INTRODUCTION 


1 


. i.  Statement  of  the  problem 


Find  the  solutions  (x), . . . ,yR(x)  of  an  ordinary  system  of  first  - 
order  differential  equations 

y\’  f‘1(*»y1»...»yn) 

(M)  . 

. *„> 

which  at  xQ  assume  n  specified  inicial  values 

(i.2)  *!(*„)-  yi0  (w . «) 

Here,  f^(x,y.j,...,y  )  are  given  functions  of  the  variables  x,y.j,...,yn. 
Defining  the  vectors 


V1 

!c  i 

. 

. 

f  « 

• 

• 

• 

• 

y 

I  f 

‘  n  / 

1  n 

we  write  (1.1)  as 

(1.3)  y'«f(x,y) 

’"e  shall  keep  to  this  way  of  writing  in  what  follows.  Speaking,  for 
-.xaraple,  of  "the  solution  y(x)"  we  mean  that  this  ist  the  solution 
vector  I  y1 (x) ' 

y(*)-|  • 

\/n<x)/ 

i.e.,  "the  solutions  y1 (x) , . . . ,yn(x)"  ,  etc. 

"'/hen  stated  as  above,  our  problem  is  already  quite  general  because  any 
explicit  higher-order  differential  equation  or  system  can  be  re-writ¬ 
ten  as  a  first-order  system.  This  requires  merely  that  all  derivatives 
exoept  the  highest  be  replaced  by  new  auxiliary  functions  (cf.  Erwe  /ll/, 
P.  27). 


1.2.  Step-by-Step  Continuation  of  5olutions 


'll  methods  discussed  in  the  following  give  reliable  approximations 


2  1.3. 

£(x)  only  in  the  near  neighborhood  of  the  initial  value  xQ.  Large 
values  of  |x-xQ|  may  soon  lead  to  poor  results*  What  one  oan  do  is 
choose  a  certain  "step"  h^  and  trace  the  approximation  only  to  the 
point  x1«xQ+h^ .This  approximation  y(xQ+h)  will  then  serve  as  the 
initial  value  of  a  new  step  from  x1  to  Xg^Xj+l^t  and  so  forth. 

/.part  from  the  specified  initial  valuef  such  "one-step-methods"  do 
not  use  any  other  information  on  the  previous  shape  of  the  solution. 
Therefore,  we  need  no  longer  bother  to  number  the  steps  but  may  oall 
any  initial  point  xQ,yo.  The  problem  left  for  the  following  chapters 
Ir;  now  to  construct  an  approximation  y(x)  at  the  point  x-xQ+h  from 
t,iven  initial  value®  xq,j  and  a  given  step  size  h  with  a  sensible 
"olume  of  calculation  in  such  a  '.ray  that  this  approximation  is  as 
close  as  possible  to  the  unknown  solution* 

I*3»  Error 


The  size  h  of  the  steps  depends  above  all  on  the  desired  accuracy. 

mailer  steps  give  better  aocuraoy  (not  considering  rounding  errors) 
but  require  more  work.  To  make  a  sensible  choice  of  the  step  size 
we  must  therefore  have  a  rough  idea  of  the  "local"  error  committed 
during  u  step  of  integration.^  shall  discuss  this  when  dealing  with 
the  different  methods  individually.  However,  the  total  error  committed 
after  several  stops  is  still  undetermined.  This  error  may  soon  become 
much  greater  than  would  be  expected  because  of  the  insignificant 
local  errors.  The  decisive  factor  is  whether  the  solutions  next  to  y(x) 
apjDroach  y(x)  or  depart  from  it  as  x  increases,  i.o.  whether  the 
solution  is  stable  or  unstable.  More  information  about  this  can  be  got 
(rom  the  so-called  transfer  matrix. In  Section  I 1.6  we  will  see  how 
t: o  calculate  it. 

!n  the  case  of  n«1,  i.o.,  one  differential  equation,  only  half  of  all 
oases  give  unstable  solutions.  In  systems  of  differential  equations 
(nonce, also  in  differential  equations  of  higher  order) ,hoacvor,  there 
is  nearly  always  at  least  one  unstable  component.  Therefore,  accuraoy 
must  be  high  should  the  solution  be  continued  over  a  domain  of 
considerable  extent.  Here  are  two  examples  that  involve  some  trouble: 


1.3- 


3 


y*'»10y  +11yv  y(0)»1,  y'(0)»-1 

y(3)-? 

(For  greater  detail  see  Collatz  / 7 / ,  p.  49) > 

y"+(l-x2)y«0,  y(o)*1 ,  y'(0)»0,  y(lOO)-?. 

In  the  last  example ,  accuracy  would  have  to  be  5000  places  if  some¬ 
thing  should  be  obtained  for  x»100. 
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> 

J 

I 

y 

Chapter  -II 


Power  Series 


by  0.  Wanner 


Abstract; 

Solving  ordinary  differential  equations  by  power  series  expansions 
has  again  become  rather  popular  lately,  on  the  one  hand  because  the 
coefficients  of  the  solutions  can  be  calculated  by  computer  through 
recursion  formulas,  and  on  the  other  hand  because  estimation  of  error 
is  relatively  simple. 


II.  1 . 
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II. 1.  Solution  by  Power  Series  Expansion 


Power  series  of  the  solutions  y(x)  will  henceforth  be  written  in  the 
way  adopted  by  ’•••'.  Groebner.  This  will  turn  out  to  be  very  useful, 
especially. in  later  chapters. 

Let  F(x,y)  be  an  analytic,  function  of  the  variables 
Inserting  solutions  y,(x) , . . . ,y  (x)  in  the  place  of  y, , . . . ,y„  -ee  find 
a  function  that  depends  on  x  only.  By  the  chain  rule.,. .its -derivative 
with  respect  to  this  variable  is 


(1 


,  f  dy,  dy  /I 

.T)  trP(x,y(x))  «  |  VFyi  "  •••  +Fyn  dTjXf 


y(x) 

Here,  the  bi’acket  symbol  £  Jx  y(x)  means. that  the  variables  x  and 
y  must  be  replaced  by  the  functions  x  and  y(x)  after  the  partial 
differentiations  have  been  performed.  Jrom  now  on  we  shall  keep  to 
this  way  of  writing,  i.e.,  every  time  some  kind  of  expression  stands 
after  such  brackets  it  must  be  inserted  for  the  variables  x  and  y. 
Since  the  functions  y{x),  which  we  have  inserted  in  Eq.  (l.l),  are 
supposed  to  be  the  solutions  of  (l.I.l)  or  (l.l. 3)  we  haver 

dyi 

=  fJ^yU))  -  I  f. 


k  •-  [if  *b<*-y>lrr-"+fn(,t’3r)ff‘l 

-  n  J 

[lPJx,y(x) 


(1.?) 


x,y(x) 


where  we  have  defined  the  linear  differential  operator 


(1.3) 


L  =  —  +f  (x,y)-—  +...+f  (x,y)-— 
3x  1v  ,'//oy.  n.  ,J/o y 

j  1 


for  brevity. 

dy  iteration  of  ( 1 . 3 )  we  find. for  the  higher  derivatives 

(1t?.)  $!L  P(x,y(x) )  =  f  DUF  1 

V.  J 


x,y(x) 


where  Ir  F  means  that  the  differential  operator  has  to  be  applied 
u  times  to  F. 

rhus,  the  ao^er  series  of  the  functions  f’(x,y(x))  at  the  point  x^ 
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II. 2. 


can  be  written  as 
(1.4)  F(x,y(x)) 


owing  to  y(xQ)-y0 
Setting  F(x,y)«yi 

(1.5)  yi(x) 


(1.1.2). 

we  obtain  the  scries  for  the  solutions  proper 
f-(x-xj"  r  1 

-Z_nrr-  {*\\  . *>• 

u-o  *-  Jx0,y0 


similar  expressions  occur  also  in  tho  theory  of  transformation  groups. 
Therefore,  such  series,  especially  the  ones  derived  in  the  following 
chapters,  are  also  named  Lie-series . 


II. 2.  Recursive  Calculation  of  the  Coefficients 


'7o  shall  now  discuss  the  recursive  calculation  of  the  power  series 
coefficients  as  lately  adopted  by  Gibbons  / 1 Q/,  R.'1).  Moore  / 36/  and 
many  other  authors.  It  has  become  very  important  through  tho  use  of 
electronic  computers. 

We  assume  that  the  functions  f  ^ (x,y)  have  been  composed  of  the  variab¬ 
les  x  and  y1#...,yn  by  finite  sequences  of  elementary  operations. 

We  note  all  intermediate  results.  Each  of  these  intermediate  results 
follows  from  one  or  two  of  the  preceding  values  (one-place  and/or 
binary  operations)  or  from  x,y1,...,yn  or  from  a  constant  c.  Suppose 
p(x,y1 , . . . ,yn) ,  q(x,y1  ,...,yn)  and  r(x,y1 , . . . ,yR)  are  three  (or  two) 
operands  that  are  interrelated  through  an  arithmetic  operation 

<2>l)  r(x-y1 . yn>‘  »(*-y1 . yn)*«<x'y1 . y„> 

or  some  sort  of  elementary  functions  g 

(2.2)  r(x,y1 . yft)-  g(p(x,y1 , . . .  ,yn) ) 

Then  we  introduce  tho  following  notation 


II. 2. 
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(2.3) 


» 


» 


. 


Hence,  these  quantities  are  functions  of  the  variables  x,y^,...,yn  • 

For  the  functions  x,y1 , . . . ,yn, f ^ , . . * ,  fn,c  which  are  special  cases  of 
such  operands,  we  shall  also  use  the  corresponding  symbols  X^,  Y^,... 

•••»Ynp»  P1p»  Pnu»  Cu  * 

In  what  follows  we  tabulate  formulas  which  permit  us  to  calculate  R^ 
for  (2.1)  or  (2*2)  from  the  coefficients 


Pp  *  Vi  »  Pp-2  Po 

»  Vi  ’  %-2  »••*»  %  (°nlY  for 

Vi  *  RU-2  *••••  Ro 


”’):r  r=p  +  q 
dffcrence:  r«  p-q 

roduct  :  r-p.q 

•'uotiontJ  r«  p/q 
oxp s  r=exp  p 

lo g‘.  r»  log  p 

"quare  root:  r«  p 

Q 

Constant  power.  r»  p 


R  «P  +Q 
•I  'I  u 

R  -P  -0 

UUP 


RU*  2_Pp°p-p  (p-0,1,2,...) 


u-1 


V  (V^I'rpS-p)/qo 

P  *0 

v  pSl!(“-p)Rppll-p  •  p0 


0*0 


P  4-P 


*  3  ?p,-  “y  (u-p)P  R  ]/P«  »  R  -log  P 

II  1.  u  IJ  Vi  P  u -p*  0  9  0  0  c 

-  f\  R  }  ,  R  -  V  P~ 

M  2Rq  l  u  p  u-p  J  *  0  o 


> _ (o<.-(o+i  )p)i*p^_p|'55-  .  VP 


p*b 


>C 

0 


( ° )  4) 


1  For  the  case  PQ»o  and  c  a  positive  integer,  G.Margreiter  has 
'.crived  special  formulas,  cf.  /53/  » 
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sin  and  cos:  q=  sin  p 


II. 2. 


I  -  -S  (u-p)R  P 
V  'i  P  U-P 


1  ^ 

r-  cos  p  Rp-  -  7>ao(ti“p}yY-p 

Proofs  can  be  found  in  Moore  /  }€>/  and  Wanner  /  5  */ • 


Q  -  sin  V 
o  o 

R  «  cos  P 
o  o 


'When  all  operations  that  give  the  functions  f^(x,y,j,...,y  )  from  x 

and  y1,...,yn  arc  replaced  by  the  corresponding  recursion  formulas, 

those  will  give  the  values  F.  from *X  and  Y.  ,  ...,Y  if  all 

1U  4  i  f  4 

derivatives  of  lower  order  are  still  present.  Because  of 
(2.4)  J)yA-  ft  ,  DU+Vi~ 

these  quantities  are  equal  to  F^=  (u+l)Y^  ^  .  Hence, 


(2.5) 


Y  a  — T.1 
li,u+1  p+1  r  i|u 


(  P*o, 1 , 2, . . . ) 

(  i*=1 ,2, . . .  )  * 


The  procedure  can  nov;  bo  repeated  with  the  quantities  Y^  •  It 
will  lead  to  a  recusive  computation  of  the  Y^  .  Recursion  begins 
with  the  values  Yio»yi  (initial  values)  using  the  formulas 

(2.6)  XQ-x,  =1 *  X2*X3",,."0 
and,  for  a  constant  c, 

(2.7)  C0“c»  0 

Then,  it  proceeds  according  to  the  pattern 


X  ,Y.  , . . . ,  Y  ...  •  f  ,  Q  ^R  ...  -*F F 

0 7  1 0  no  0 ’  0  0  1 0  no 

_ 4. _ t _ > 


Yn»M,,Ynr  v** 


nl 


Y12»*,#,Yn2  * 


•••R2’^2  ^2*" 


1 2 1 


< - ~r 

Ili».»»»Yz  etc. 


n2 


1 3»#,,»4n3 

Some  of  the  authors  that  have  worked  with  one  of  these  (or  similar) 
recursion  formulas  are  Steffensen  /46/,  Miller-Hurst  /35/,  S.Rabo  /39 /, 
V.Oautschi  /16 /»  E.Fehlberg  /14/,  I. iennig  /  /,  Doprit-Zahar  /  9/, 
Leavitt  / 32/ >  Richtmycr  /41/. 


II. 3 .  Estimation  of  Error 


Estimating;  the  error  of  a  oeries  that  hue  boon  out  off  (e.g. ,  after 
the  m-th  term)  is  thus  indispensable  for  a  sensible  ohoioe  and  control 
of  the  stop  size.  One  possibility  is  to  bound  the  error  by  means  of 
majorant  series  (e.g.,  W.  Groebncr  / 2l/,  /  22/)  • 

However,  with  Buffing’s  differential  equation  as  an  example,  G.  Maefi 
/  34/  has  shown  that  using  Lagrange’s  remainder  of  Taylor's  series 
gives  an  error  limit  which  is  by  3-4  powers  of  ten  more  accurate  than 
in  the  case  of  the  majorant  technique. 

If  ->.11  occuring  derivatives  exist  and  are  continuous „  then  we  have, 
according  to  Lagrange, 

yA(x)  - 


with 


im 


m+1 
(m+1 ) ! 


(x-xo) 


,n+1 


iv. 


m+1 


y,M 


•  Vf‘* 


where ,  owing  to  (1.2') 
(3.1)  R 


(x-*J 


m+1 


im 


o' 
(n+1 ) ! 


$  »y(f) 


,  x0^Jx 


For  a  precise  estimation  of  the  error  one  has  to  know  a  domain  B  wJ.ic. 
is  known  to  contain  the  solution  y(|).  The  functions  Da+1yi  can  then 
be  estimated  in  this  domain  (Fig.  l) 


y  t 


i 


\ 


i 


(Pig.  1) 


i  i 


i/'aeu  /$4/  demonstrates  this  by  Buffing's  differential  equation. 
fl.E.  Moore  /36/  solves  this  problem  generally  and  automatically  by 
m^nns  of  interval  arithmetics- 


12 


II. 4- 


..noving  an  approximate)  «rror  is  sufficient  for  a  sensible  control  of 
the  step  size.  Here,  one  may  put  up  with,  3ay,  the  value  of  D'+1y.» 
at  the  point  xQ,yo  (this  would  be  the  first  term  neglected),  or  rather-, 
one  chooses  the  larger  one  of  the  values  ut  the  points  x0tY0  and 
xo+i.  ,y(x  +h)  (starting  point  for  the  subsequent  3trp).  Both  numbers 
arc  easy  to  computes  it  is  sufficient  to  run  the  iteration  for  calcu¬ 
lating  the  Taylor  coefficients  for  this  and  the  next  step  through 
anothor  loop. 

Vo  shall  obtain  the  formula  (  3  . 1 )  for  the  remainder  as  a  special  case¬ 
in  the  next  chapter. 


TI-4.  Transfer  Matrices 


Lot  yA(x)  bo  solutions  of  the  differential  equation  (i.I.l)  for  the 
initial  values  ykQ.  The  matrix 

/  9yt(x) 


(4.1) 


H(x) 


*  (Hikw)  ■ 


3y 


ko 


.’hich  consists  of  the  derivatives  of  the  solution  y^(x)  with  rospeot  tf 
'ho  k--th  initial  value  yko,  is  then  called  the  transfer  matrix  per¬ 
taining  to  y(x). 

In  other  words;  The  transfer  matrix  describes,  in  first  approximation, 
the  variation  of  the  solutions  y^  at  the  point  x  if  the  initial  values 
y .  arc  changed.  When  we  change  the  initial  values  y.  , . . . ,y _ by 

i'^O  J  Q 


no 


,  ...,cno»  the  solutions  at  the  point  x  will  in  first  approximation 


-1o 
change  for 


e  Ax) 


*yA*) 


9y 


1 0 


C4  + 
lO 


ay^x) 


Thus , 
(4.2) 


/  ei(x> 


1 

VnMl 


I  Sy^x) 
dy1o~ 


\  6yn(x) 
\  ?y1o 


9y, 


no 


no 


Sy/x) 

3y. 


no 


fyri(k) 


jy 


no 


M 


\tnol 


or,  in  vectorial  form, 
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(4.2') 


x 


For  example,  errors  committed  somewhere  in  the  numerical  integration 
can  be  mapped  forward  or  backward  to  any  fixed  point  by  moans  of  the 
transfer  matrices. 

Hence,  these  also  describe  hov'  an  error  committed  at  a  certain  place 
influences  the  final  result,  ^e  shell  consult  the  transfer  matrices 
also  for  an  "optimum"  step  size  control  which  takes  stability  and  the 
total  final  error  into  acoount  (Chapter  VI). 

i’ho  transfer  matrices  arc  also  useful  in  boundary  value  problems  in 
which  some  of  the  initial  values  are  missing  and  have  been  repluood  by 
conditions  at  other  parametric  points.  Hero,  the  missing  initial  values 
must  first  be  guessed  and  then  bo  improved  by  mcanB  of  the  transfer 
matrices,  after  tho  relevant  solutions  have  been  found  (Y/anner  /3V)» 


f 
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A  remainder  formula  for  the  error  after  the  m-th  term,  which  is  analo¬ 
gous  to  (3*1),  is 

z0ilix  ■ 

In  the  case  of  a  step-by-step  integration  of  the  differential  equation 
with  tho  intervals  xq<x1<  . . .  < x^  ,  Eq.  gives  only  the  looal 

transfor  matrices 


3yi^xi^ 


Owing  to  tho  chain  rule  (for  functions  of  several  variables),  these 
matrices  can  be  multiplied  with  each  other  to  give  the  total  transfor 
matrix 

( r  •  3 )  H(xn)  -  C(xN,xN_1)...C(x2,x1) C(x1,xQ) 

Notice  that 

(5.4) 


H(xq)  «  C(x,x)  «  E  (Identity  matrix) 


and 

(5-5) 


C(x,x* )  -  C(x' ,x) 


-1 


For  linear  systeris  of  differential  equations,  the  columns  of  the  trans¬ 
fer  matrix  coincido  with  tho  fundamental  solutions  of  the  oorrespondir.e 
homogeneous  system  (with  the  initial  values  0, . . ,  1 , . .  ,0) ,  and  the  rela¬ 
tion  (4-2)  not  only  holds  in  first  approximation  but  is  valid  exactly. 

Another  possible  way  of  calculating  the  transfer  matrix  is  to  integral^ 
tho  system 


dHik<x)  V1" 
~ 35  "  fct 


•Ji 

dy 


J  J 


x,y(x) 


Vx) 


for  every  k«1,...,n  with  the  initial  values 

Hik<xo>  •  bik 

together  with  Eq.  (i.I.l).  This  formula  is  usually  given. 


II. 6 
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ir.6.  Recursion  Formulas  for  Expressions  with  an  Additional  Operator 


Here,  we  replacotho  operator  6 f  (5.1)  generally  by  D  because 

v/o  chall  need  the  following  formulae  for  other  purposes  too. 

Suppose  D  is  another  linear  differential  operator.  In  addition  to  (2.3  • 
we  adopt  the  symbols 


(6.1) 


r  Ml 
V  p! 


for  certain  operands  p,q,r.  Again,  these  quantities  are  functions  of 
Xvy1.-**»ynt  ;tnd  the  corresponding  symbols  X  , Y.,  , . . . , Y  , . . . ',P  , 

are  again  valid  for  the  functions  x,y1 , . . . ,yn,f1 , . . .  ,fn,c. 

Also  for  these  quantities  we  obtain  rocursion  formulas  by  simply 
applying  the  operator  D  to  the  formulas  of  page 


Table 


Sum;  r*p+q 

R  «P  +Q 

u  U  u 

Difference;  r«p-q 

vv% 

Product  r=p.q 

_  -  •  - 

Rp"  ^p^U-p+?p%-p 

Quotient:  r»p/q 

R  ■’{  P  -  /  (R  Q  +R  ft  )-R  Q  I/O. 

^  t  jj  /  '  p  p  p  V|l»p/  (|^o  0 

P  *  0 

oxp:  r*=exp  p 

R  =  -  >  (u-p){R  P  +R  P  },  R.RP 
u  u  L_'  p  u-p  p  p.-p-'’  0  00 

log  r»log  p 

R>|P  ("-p)(P  R  +P  R  )-P  R  ]/P 

U  N  U  (1  — -  '  r/'  P  u-p  P  u-p'  0  \1>  0 

R  nP  /Pn 

0  0  0 

11=1 

Root;  r=  p 

Rp=  2R”  ( Rp"2i—  RpRp-pJ  P-0,1,... 

0  p-o 

Constant  power t  r*pc 

VU^I(cu-<<!+l)|,)(5pVp+RpVp>-%fo}/p 

v° 
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:  in ,  ccs  •  q=sin  p 

r=cos  p 

Si.'c  Wanner  /51/,  p«  27* 


II .  7 . 


±zl 

T\  -  -  /  (u-pJIt?  p  +RP  ?  .  ,  “ 

u  4  p  p-p  p  4-P  J  o  o  c 


•t-1 

R  =  —  )  (u-p)(q  P  +Q  P  \  >  R  --Q  ? 

|i  4  /  _  vr  r,C  p  |i-p  4«p  )  7  o  00 


First  of  nil,  all  expression  must  exist  should  R^  be  calculated, 
.pp  lying  1)  to  (2.5)  we  obtain 


<  r. 


2) 


*1,4+1 


which  enables  us  to  employ  recursion.  As  wo  can  soe,  the  above  formulas 

pro  independent  of  the  particular  choice  of  the  operator  D.  Setting, 

d 


0  •  g  • ,  D  = 


3y„ 


we  find  the  expressions 


Y. 

14 


1 _ o_ 

Ml  ayk 


which  arc  needed  in  ( 5 . 1 )  •  In  thio  case,  recursion  starts  with  the 
initial  values 


For  the  independent  variable  x  ve  have  XQ  =  X1  »=...*  0 

and  for  a  constant  c  C  «  C 0  . 

0  1 


Subroutines,  which  calculate  these  formulas  are  given  in  Knapp -Wanner 
/30/  or  Wanner  /5I  /. 


11.7^  Recurs ion  Formulas  for  Other  Operations 


The  class  of  operations  th '.t  are  allow?!  for  the  formation  of  the 
functions  f^(x,y)  will  be  consi1  rably  expanded  in  this  section.  We 
shall  3oe  that  every  function  satisfying  a  differential  equation  that 

cart  already  be  processed  is  allowed  her". 

First >  vu.  show  by  way  of  a  few  examples  how  recursion  formulas  can  be 
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got  for  many  functions  by  introducing  auxiliary  expressions  i 
rgqrctan  p  : 


Here  > 


Dr 


Up' 


V7 o  sot  1+p  *q  ,  whence  qDr*Dp 
Wc-  find 


q  ■  y~ P  P  ,  Q  «1+P  ‘ 
4  fso  p  M"p  0  0 


(7.1) 


R 


-{p  -  -  V_(u-p)0.  R  ]/0  , 

U  V  U  U  *■ — r  0  U-pJ'  0  * 

0*1 


R  =arctan  P_  . 
o  o 


r=tan  p  When  sin  p  and  cos  p  occur  simultaneously,  the  best  way  is 


to  write  r  =  an^  to  uso  the  formulas  of  page 

cos  p  r  ® 


sin  p 
cos  p 

7hon  sin  p  or  cos  p  docs  not  occur,  it  is  preferable  to  U3e  the  formulas 

2 


JLzl 

■r2_J 


Vr^Vn-i  ’  V1+Ro 


(7.2) 


R  »  —  )  (ji-p)Q  P  ,  R  =  tan  P 
4  4  P  u-p  ’  o  o 


which  have  been  obtained  by  reversing  the  formulas  (7*1 )• 


r=arc3in  p  : 

i  A  p  p 

Wo  sot  q=V(l-p  )  ,  whence  qDr*=Dp  .  Owing  to  q  *1-p  wo  obtain 


V> 


(7.3) 


V^^pVp^SVp]  • 

V  l  V  jS'-PlVu-pl^o  ’  Ro'arc8in  Pt 

P=1 


For  r=arcco3  p  ,  all  formulas  remain  the  same,  except  for  R0«arccos  T>r 

For  the  corresponding  hyperbolic  functions,  only  a  few  3igns  have  to 
bo  changed  in  the  formulas  on  page  16.  Of  course,  also  for  all  these 
formulas  there  are  also  the  corresponding  recursions  with  the  additio¬ 
nal  operator  D. 


Consider  the  general  case  that  (x) , . . . ,um(x)  are  solutions  of  the 
differential  equations 

u{(x)  -  gi(x,u(x))  . 

The  only  assumption  we  make  is  that  the  funotions  g^  are  m  ide  up  only 
of  the  operations  dealt  with  so  far.  Then  we  can  give  recursion  formu¬ 
las  also  for  those  functions.  This  step  can  be  repoated  over  and  ovo  • 
and  loads  to  a  successive  extension  of  the  recursion  formulas  to  more 
and  more  functions  of  analysis. 

Let 

r*  -  u^p)  (i»1,...,m) 


Since  the  functions  gi  are  made  up  of  operations  whose  recursions  arc 
known,  wo  can  calculate  the  expressions 


* 

i 


,U-1 


D^gjCPtuCp)) 

(u-i)l 


for  the  operand  p(x,y)  from  the  coefficients  up  to  P^,  ^  fJ_1 
Because  of  Dr^-u^(p)Dp  =  g^(?»u(p))Dp  we  have 

(7-4)  R .  -  -)  (u-p)G*  P  ,  R.  -u.(P  )  , 

v  '  iu  |i4 — _v  ip  u-p  ’  io  iv  o' 

p= 0 


the  sought  recursion  formula. 

Finally,  we  consider  the  important  equation 

(7.5)  *,(x)u"  +  a2(x)u‘  +  n.j(x)u  »  0 


v/hich  with  the  usual  substitutions  u«u<1  ,  u’  =  u2  becomes 

ui  *= 


u* 

? 


ai  Wi? 


Let  p(x,y)  be  an  arbitrary  operand  and  lot 
r1  =  u1(p)-u(n)  ,  r2=u2(;:)«=u'  (p) 


We  put  ^jc(p)  =  ^(x,y)  and  assume  that  the  coefficients 
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4kp 


DPak(p) 

P* 


can  be  calculated  by  moans  of  the  existing  recursion  formulas  from 

the  values  P  ,  ...,P  . 

o’  ’  p 

Moroover,  v/e  use  the  notation 

*  *  4 

«1(p)u1(p)  +  a?(p)u2(p)  -  +  ajr2  -  q  ,  -j  -  a  ; 

a3 

then  we  find  the  recursion  formulas 

QM-l“^l-A1pR1.M-P-1+A2pR2.U-P-1-i  ! 


p-o 


n-l"  tS-l"  ^SpA3.4-p-l-i^A3.< 


(7.6) 


R  -  -  /  (ji-p)R  P 
1H  4  ?P  U-P 

R  =  ”  /  (u-p)S  P 

2u  |iZ — v'  o  u-p 
p*o 


J 


(p=1,2,...) 


R1,o-  “<Po>  -  *2,0"  “’<Po>  • 


I'nesc  recursions  are  valid  for  all  functions  that  satisfy  a  differen¬ 
tial  equation  of  the  form  (7*9) >  that  is,  for  example,  all  kinds  of 
Bessel  functions.  Mathicu  functions,  Weber  functions.  Chebyshcv-. 
Legendre- ,  Hormite- .  Laguorre-.  or  Jacobi  polynomials .  etc. 
for  Bessel  functions  of  the  first  kind,  e-g. ,  wo  have 


-  P 


2 


l2. 


l1 


_  2  2 
p  -n 


and 


A,  -  7  P  P  , 

?M  4— r  P  U-P  * 

p  =  0 


74 


*  *  0 

A.  =  A,  -n*6  ^ 
1  u  3u  uo 


Finally,  we  should  like  to  mention  that  for  orthogonal  polynomials, 
in  particular  for  higher  n  (low  n  are  uninteresting),  the  above  for¬ 
mulas  require  much  less  work  than  using  the  "generating  functions" 
as  suggested  by  Leavitt  2/. 


LIE-Series 


by  G.  Wanner 


This  chapter  discusses  the  numerical  evaluation  of  W.  Groebner's 
Lie  series  perturbation  formula,  on  v/hich  an  efficient  numerical 
method  v/ith  satisfactory  error  estimation  is  based. 


III. 2. 
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III . 1.  Groobner's  Perturbation  Formula 


Groobner's  perturbation  formula  states  ho v/  one  has  to  correct  an 
arbitrary  given  approximate  solution  y(x)  ($.j  (x) , . . .  ,yn(x) )  in  order 
to  find  the  solution  y(x).  This  formula  is  a  generalization  of  Taylor'? 
series  (II.  1. 5)  which  can  be  obtained  from  it  when  the  operators  are 
choaon  in  a  special  way. 

A  3y3tem  of  differential  equations  must  be  known  for  the  approximate 
solution  y(x): 

(1.1)  y'  -  f(x,y) 


and  the  approximate  solution  must  assume  the  samo  initial  values 
(1*2)  y(xQ)  *  yc 

:’rf;  introduce  the  operator 


(1.3)  *  ( f  1  ( x , y )  f^x.y))—-  +  ...  +  (fn(x,y)-?n(x,y))*^- 

1  Jn 

which  accounts  for  the  difference  between  th>-  two  differential  equation 

Ilono  G  y 


D  -  D1  +  D2 


(1.4) 

As  ’./o  shall  see  in  the  next  chapter,  here  we  have  the  formula 

d.5)  y<*>  *  9M  +  II  |  [v"0|,Jq)  4f 


('7.  Groebner)  for  the  sought  solution  y(x). 


III. 2.  Knapp's  Remainder  Formula 

According  to  Knupp  / 26/ ,  the  remainder  of  the  series  (l  5)  after,  say, 
s  terms  it, 

(2.1)  y(x)  =  y(jc)  +  R,„(x) 

■J 

with 
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S  X 


2-1')  H*)  -  y(x)  +  L_  £  C1>2D“yJ f ^ 


(2.v)  *.(«)  - 1  ^|[»',*V],,y(?)-lj>S+1y]?,5(f)ld!  ' 

o 

We  shall  prove  this  formula  in  the  next  chapter.  Knapp  /  26  /  has 
derived  these  formulas  assuming  that  f^,?^€C3  . 

Another  formula  for  the  error  can  be  obtained  by  increasing  s  in  (2. 
by  1  and  adding  the  last  term  ot=s  +  1  of  ( 2 . 1  * )  s 

,  (x. c's+1 

(2.2)  R  (x)  '  ^ 


^l&^+VJ4,rt,rk*+V]|>W|)]d, 


Finally,  a  mean  value  theorem  of  integral  calculus  gives 

.3  +  1 


(x-xj 

(?.T>)  R  (x)  . 


(a+1 ) 


—  [[,3+1Sr]  ^  ,y(| ,  )-[?’*  ^ 


*04  ?14* 


and 

(2.4)  Rs(x) 


(x-xJ 


\s+2 

f-— 1&i>3+v]5  (§  P(f 


(a+2)l 


V  MX 


III.?.  Special  Cose-  Power  Series 


The  power  series  expansion  of  the  solution  y(x)  is  a  special  case  of 
the  Lie  series  (2.1),  if  the  origin  .1  differential  equation  is 
autonomous,  i.  e. ,  if  the  function  f(x.,y)  do  not  dopend  on  x  and  vo 


have  D  =  t—  +  f(y)|~  • 

ox  v,7/dy 


To  show  thi3  we  put  D.  = 


d_ 

Ox 


=  f(y)§^  »  y(x)*y0  ,  thus 

since  hero  nl3o  Day  do  not  denend  on 


25, 


How  the  integrations  are  readily  carried  out  giving  with  a+1«B  and  the 
remainder  (2.4) 

7f~\°  yV<»*>  '  -» 


,00 .  ft 

B-o 


o  (8+2) 

tho  fonaula  (ll»3*l)  v.’'-»ich  has  been  used  before 
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III. 4.  Choice  of  Approximate  Solution-: 

Tho  approximate  solutions  y.,(xj , . . .  ,yf  (x)  cur.  he  chos  en  freely.  They 
only  have  to  satisfy  the  initial  conditions  (1.2),  end  a  syotom  of 
differential  equations  must  be  knowr.  for  them.  The  bettor  tho  choic^ 
of  the  approximate  solutions,  th  ore  efficient  is  the  method. 

It  io  oxpidiont  to  use  the  first  terms  of  tho  power  series  expansion 
(II. 1.5)  for  an  .approxir.*  tion  to  scare  "ith: 

(4.1)  V*>  ’  tl  “IT*“ 'P"7! j*0./0  ' 


u»o 


~o’yo 


(of.  (II..7. 3))i  of  course,  n  may  *).  o  uepend  on  i.  A  corresponding 
system  of  differential  equation?  con  bo  found  by  simply  dif  forontiuti.i 
Kq.  (4  1)  (the  quantities  |  Y  “I  are  constants) 


(4*1 ' )  ? 


(cf.  (II. 2. 5)),  whore  tho  functions  depend  on  x  only.  The  formulas 
(4.1)  r.nd  (4.1*)  are  used  in  she  general  nrorran  CR0*!BII5!I,  reproduced 
in  Knnpp-Wanner  /30/  or  banner  /51/. 

One  may  also  retain  parts  of  .h •  ■  original  system,  e.g. ,  in  eouationn 
of  the  kind 

yi  -  y2 

y2  *  y3 


mm 
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and  replace  only  the  last  equation  by  a  polynomial  in  x.  In  this  case, 
however,  the  degrees  m  in  Eq.  (4<l)  must  decrease  by  unity  each  ns  i 
increases.  This  reduces  the  operator  D?  to  a  simpler  form  since  then 
it  consists  merely  of  a  single  term. 

When  the  functions  Dm+1y^  are  bounded  in  a  region  B  of  (x,y)-space, 
then  we  have  from  ( II . 3 • 1 ) 

|x-x  |m+1 

(4.2)  ly.M-y.M  I  4  0  i— - 

1  1  (m+1 ) ! 

or 

(4*2')  *  0((x-xQ)ra+1)  . 

In  this  case  we  say  that  y^(x)  is  of  the  order  m  (or  of  the  error  ordc: 
m+1 ) . 

Of  course  there  are  examples  for  whioh  n  choice  other  than  (4*1 )  is 
more  convenient,  e.g. ,  the  equation 

y'  -  Vx  +  V"y  ,  y(0)  »  0  . 

Here,  the  first  Taylor  term  vanishes,  the  second  is  infinite.  However, 
choosing  the  approximate  solution 


wo  obtain  from  (2.1')  with  s»0 

y  -  f  . 

Compared  with  other  methods,  this  in  a  very  good  approximation  /4?/. 
For  small  x  values  it3  accuracy  is  sufficient  and  the  singular  point 
x*0  can  be  avoided.  More  terms  of  Fq.  (2.1')  are  not  allowed,  because 

a  4 

only  f€C  ,  whereas  f$C  . 

III. Order  of  the  i-Iotuod 

Definition- 

A  method  io  of  tho  ordor  p,  if  the  solutions  $(x)  obtained  through 
it  are  of  the  order  p,  i.g.,  if  for  every  solution  y(x),  whose  Taylor 
series  exists  far  enough, 
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yt(x)  -  ^(x)  *  0((x-xQ)p+1y 


The  Taylor  series  of  the  tv/o  solutions  \i ill  then  agree  up  to  et  leaf -i 
the  p-th  term.  •*  • 

We  shall  prove  now  that  the  method  dofi~od  by  Eq.  (2.1')  is  of  the 
order  m+3+1 ,  if  the  otarting  solution  y(x)  is  of  the  order  m: 

Theorem :  If  the  functions  Daf^  oatisf-  ;  Lipschitz  condition 

(’•’)  ID>%]  .-[»%]  Lly*  -  y**| 

v  ir  ^  -*  v  v  ](=  ! 


x,y 


x,y 


4  M 


K  n|x-x 

O 


'  m-  :+2 


in  a  region  _B,  and  if  the  starting  ac 7 ution  is  of  the  order  m 

I  I  n+ 1 

(5.2)  |y±(x)  -  y4(x)  |  &  M  — - -  > 

1  1  (a+1 ) ! 

than  the  relation 

(5-3)  ly^x)  -  ^(x)! 

1  (m+oJ  )' 

holds  true  for  the  solution  y(x)  thoi  f  ~Mows  from  (2.1 ')» 

Proof.  D0+1yi*PDf^  »  hence  3ubstitutd  <■*  5tl)  and  (5.2)  in  (?.1" 
get 

I yi(x)  -  yt(x)  |  -  |  RiJ3(x)  |  - 

i|  “iW-vj,, »J, «  ,UU 


•l|  V  lac  ,  I 


iio"  the  statement  follows  by  means  of  th<-  well-known  integral  foyr’ula 
iU  x  /„  „  \a  /.  ..  .*• 

V  A**  A  _  j 

(5*4) 


(x-xo)' 


f  (x-l)a  (?  -X  V  - 

i  —  V~  rr  di  > 


u! 


This  theorem  is  a  special  ca^e  of  1  *  u  al  theorem  stated  in  th< 
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next  chapter. 


Thus,  the  order  of  the  method  increases  by  1  with  each  additional 
integral.  It  may  also  increase  by  more  than  1,  as  for  example  in  the 
following  case; 


,  10C>  100 

y»  -  x  +y  , 


L-  +  (x100+y100)- 
dx  KX  y  '  dy 


y(o)  -  o 

.100  3 


\o  _  o  iuu  o  _  100  9 


101 

Hr)  - 


where  Eq.  (2.1')  with  s«1  gives 

A  }  .10100  5 

y(x)  -  y(*)  *  J  '|Q— Too  d!  * 


100 


20099 


(101) 


199 


d|  + 


x101  X10101  x 20101 

101  10101 (101 ) 100  20101. 201. ( 101 ) 199 
this  contains  already  more  than  30  000  Taylor  terms. 


III. 6.  Numerical  Evaluation >  Quadrature  Formulas 


To  evaluate  Eq.  (2.1')  numerically  wo  must  calculate  the  occuring 
integrals  in  a  proper  way.  Tho  following  lemma  is  quite  useful  for 
this  purpose. 

Lemma:  If  the  starting  solution  y(x)  is  of  the  order  m 
yA(x)  -  y*(x)  -  o((x-x0)::l+1) 
and  if  the  fi(x,y)  satisfy  a  Lipschitz  condition,  then 


(f.l)  fk(x,$(x))  -  fk(x,y(x))  *  0((x-xQ)m)  . 


Proof:  According  to  tho  first  assumption  we  have 

/}(*)  -  y}(*>  -  o ((x-*o)“) 


*) 


because 


yi(x)-yi(x) €  cm+1 
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and  owing  to  the  Lipsohitz  condition  we  have 

f\(x,y(x))  -  fi(x,y(x))  -  0((x-xo)m+1) 

Hence , 

fi(x,y(5:))  -  i\(x,y (?:))  =-• 

4s  t  A 


«  f^x,^*))  -  (x,y(  k))  +  r.(x,y(x))  -  fi(x,y(x)) 

-  0((*-xo)m1) 

*  0((x-xJm) 


+  yj(x) 


-  y{(*) 


Nov/  ufc  have  to  calculate  the  following  integrals 

X  X 

<«•*>  i  [v\],  „,„/•<  '  J  «<!>*!  • 

*o  '  *o 

We  choose  the  step  size  h  and  set,  as  usual  according  to  Gauss, 
x  +h  .  k 

0  ■  i 

(*•3)  J  c 

xo  Js1  ’ 

v/horo  the  a.  (Oia4il)  determine  cho  given  basic  points  at  which  the 

•1  ij 

vnluoa  of  *(?>■  which  are  then  summed  up  with  the  weights,  must  be 
calculated.  The  root  of  this  section  will  now  be  dedicated  to  deter¬ 
mining  tho  coefficients  a.,  and  c.  ap  expediently  as  possible. 

.10 

First,  v/e  find  from  the  lemma  that  the  function  g(f)  contains  the 
factor  (|  -xe)m,  for  r;<=  have  (of.  (l.*>)) 

hjr.co,  owing  U>  (6.1) ;  ir  ,i  factor  occurring  in  the  braces. 


Thus , 


rx%)  - 


and  (6.5)  attains  t 

ho  form 

xof 

(a,)  j'5,,1 

r*r(f  )dj 

. 

-  V+1  C 

X 

0 

f\r  1 

Tv'u  transformation. 

*  -x  -r.t 

give . 

.a . 

0  o 


l  Tf".  ,1 

0  ,1 
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(6.5)  \  tQG*(t)dt  -  y  C.C*(a.) 

0  0^1  0  J 

with 

(6.6)  G  (t)  ■  G(xQ+ht) 
rind 

(6.7)  C,  -  c.a1" 

J  J  3 


Equation  ( ^ . ^ )  shows  how  tho  coefficients  and  must  be  determined 
so  that  an  order  as  high  as  possible  will  be  attained?  The  must  bo 
the  zeros  of  the  k-th  one  of  the  polynomials  whioh  in  the  interval 
(0,l)  are  orthogonal  with  tho  vveight  function  tm,  the  are  the 
corresponding  weights  (o.g.*  Natanson  /38/,  p.  436).  These  coefficient.- 
are  tabulated  with  8D,  e.g. ,  in  Krylov-Lugin-Janovich  /3V*  Stroud- 
Secrest  /49/  give  a  FORTRAN  program  for  this  (however,  for  tho  interval 
(-1,+l)).  The  coefficients  c^  can  then  be  found  by  means  of  (6.7) - 
They  can  be  calculated  explicitly  for  k=1,2; 


k-1  s  ai]  * 


m+1 

m+2 


m 


(m+l)a“ 


k»  2 .  a 


m+2±V27m+2)/(m+3) 


1 , 2  m+4 

ci  *  -  °2  37?)  1  (''“('V';-)) 

C2  ‘  ^0^2  ■  n1  !  ^2^a2*a1^  • 


III. 7«  Some  Values  of  the  Tabic  of  Coefficients 


Here  are  the  coefficients  ,  c^  of  the  quadrature  formula  (6.3)  for 
a  few  values  of  m  and  for  k*l(l)4  with  an  accuracy  of  about  23  places. 
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III. 8 


III. 6.  Effective  Formulas 


To  calculate  the  integrals  (5.2)  by  means  of  (6.3)  we  must  evaluate 
DpDay.,  at  the  point 


(s.i) 


»+ajh  *•  -  y(Vajh) 


If  v/e  want  to  do  thio  by  means  of  the  recursion  formulas  of  Sec.  II. 6. 
we  have  to  calculate  the  expressions  f~Yia  j  with  the  formulas  of 

■  otJVno 

II. 2.  first,  because  these  are  needed  for  the  general  recursion  formal 
Then,  the  formulas  of  Sec.  II. 6.  give  the  expressions  (cf.  ( II .6.1)) 


i'vM 


r—  i 
y  * 


where  (cf.  (l.3))  iteration  must  be  started  with  the  values 

[YioJ^,>2.  "  [ViJfj,^  =  fi(*j’V  “  fi(VV 


and  where  wo  have  to  put 


X  *  X.  0 

o  1 


for  the  independent  variable  x  and 


C  =  C.  =  ...  *  0 
o 


for  a  constant  c. 

How,  the  formulas  (6.2,  6.3)  assume  the  form 
x_+h 

(8.2) 


ot,h  (x  +h-l)a  _ 

1  0  ivftl 


a! 


ds  = 

a/,  \  ; 


o 
.  k 


=  h 


t  ,v(\) 

]  \j’\) 


h“+1  £•,<  wlvl. 


,.o+1 


Zv.  ft  ] 

4— r '  i.  i«j| . ,  n . 
J-1  M'\l 
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where  the  quantities 

(9.5)  ria  -  Ojd-'j)0 

can  be  prepared  at  the  very  beginning. 


III. 9*  Choosing  the  Orders  n»  a,  and  k 


For  choosing  k,  i.g.,  the  number  of  the  base  points  used  in  the 
quadrature  formula  (8.2),  it  is  important  to  consider  that  the  errorc 
of  the  quadrature  formulas  and  the  methodical  error  caused  by  breaking 
off  the  series  (2.1* )  should  bo  of  the  sane  order  of  magnitude.  Other¬ 
wise,  it  would  make  no  sense  going  through  the  trouble  of  calculating 
higher  terras  of  (2.1  * )  'while  ar.  error  tc-n  times  as  large  has  already 
been  committed  in  the  quadrature  of  the  first  (and  usually  largest) 
integral.  On  the  other  hand  it  makes  also  no  sense  to  calculate  the 
integrals  with  particular  accuracy  in  view  of  a  large  breaking-off 
error.  Wo  shall  therefore  try  to  choose  k  in  3uch  a  way  that  the  qua¬ 
drature  formula  13  of  at  least  the  same  order  a3  the  method,  but  then 
its  order  is  not  much  higher  either.  As  is  known ,  Gauss's  quadrature 
formula  with  k  base  points  is  of  the  order  2k.  By  means  of  the  lemma 
in  Sec.  III. 6  we  succeeded  to  split  off  the  factor  V  fron  the  inte¬ 
grand  (of,  (6.5)).  Therefore,  the  nrder  of  the  quadrature  formula  has 
been  raised  to  m+2k.  The  method,  on  the  other  hand,  is  of  the  order 
r.i+s+1  (cf.  Sec.  III. 5).  Equating  both  orders  we  have 

(9.1) 


Hence1,  k  should  be  about  half  as  great  as  the  number  of  integrals 
used. 

The  choice  of  m  and  s  ia  a  question  of  the  differentiability  proper¬ 
ties  of  the  difforexitial  equation  as  well  as  a  question  of  expenditure- 
With  the  quoted  recursion  formulas,  labor  is  approximately  proportio¬ 
nal  to  (m+l)m+2k(s+l)(s+2)  or,  with  (n.l ),  to  (n+1 )m+(s+1 )2(s+2)  ; 
thus,  it  increases  with  c  ranch  faster  than  with  \i.  Minimizing  this 
expression  under  the  subsidiary  eond-i  tion  of  constant  order  m+3+1 
one  finds  (s+1 )(  •;  this  "pnlics  to  the  combinations 
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s  I  0  1  j  2  ;  3  i  4 

- 1 - 4- - 1  .-~u — - - 

o  12  '  0  <  16  j  27  U2 

The  choice  of  n  and  s  is  also  a  question  of  the  desired  accuracy.  The- 

influence  of  m  and  s  on  the  results  depends  on  the  magnitude  of  the 

constants  M  and  K  of  the  theorem  in  Sec.  III. 5-  A  nethod  of  higher 

s 

order  is  always  better  than  a  method  of  lower  order,  if  the  error 
limit  is  small  enough.  U3ua  ny,  n  is  chosen  between  5  and  20,  s  between 
0  and  5-  Then,  with  the  limits  of  accuracy  chosen,  one  tries  to  fit  th 
stop  size  along  the  solution. 


III. 10.  Estimation  of  Error 


The  (methodical)  error  committed  in  one  step  may  be  estimated  by  mcanr 
of  the  theorem  of  Sec.  III. 5  (e.g. ,  by  regarding  the  difference 
y(x)-y(x)  as  the  error  in  y(x))  or  by  directly  estimating  the  remain¬ 
der  formula  (2.1").  The  latter  case  will  be  considered  here.  We  may 
replace  (2.1")  by 


(10.1) 


# 

R. 

is 


c0  s*  u  -s,y(?) 


for,  owing  to  (2.1") 


i3  the  error  in  the  solution  obtained  from  (2.1’)  when  y(x)  is  used  as 
a  starting  approximation.  But  owing  to  (5*2,  5*3)*  we  have  for  this 
error 


r**w! 

is v  ' ! 


,  |K  nlx-x^ 

*-*  )m+1  ^ - S. 

:  O' 


3+1  T 


(m+1+2(s+1 )) ! 


* 

i.e.,  its  order  is  much  higher  than  that  of  R^b  ,  therefore,  it  way  bo 
neglected . 


The  following  lemma,  which  is  about  the  order  of  the  integrand 
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i3  useful  to  an  expedient  evaluation  of  (lO.l). 

Lenraa ;  From  the  Lipschitz  condition  (5*0  for  D8f^»D8+1y^  and  from 
y^xj-y^x)  =  o((x-xo)m+1) 
it  follows  that 

zA(x)  *  O((x-xo)m+1)  (i-1,...,n). 

Proof;  From  the  theorem  in  Sec.  III. 5  follows  y^(x)-y^(x)»0( (x-xo)m+S+2) 
This  relation  and  yi(x)-yi(x)«(yi(x)-yi(x))+(yi(x)-yi(x))-0((x-xo)n+1) 
together  v7ith  the  Lipschitz  condition  give  the  statement. 


Hence,  we  have 

(10.3)  Zi(x)  -  (x-xQ)::i+1Zi{x)  . 


With  this  expression  we  approximate  the  integral  ( 1 0. 1 )  by  means  of  a 
quadrature  formula  which  uses  only  tho  point  xQ+h,  i.e.,  the  end  point 

of  tho  step  of  integration,  as  a  base  point.  The  funotion 

(x-^)S(«-xQ)n+1  is  split  off  as  a  weight  function.  Therefore,  we  put 
xQ+h  xQ+h 

(10.4)  J  (x0+h-?)sZi(|)d^  ]  (xo+h^)s(|-xo)n+1Zi(^)d|-cZi(xo+h) 

X  X 

0  0 


'vC  determine  the  weight  factor  c  in  such  a  way  that  (10.4)  is  fulfilled 
exactly  if  Z^*)  is  constant.  This  gives 


c 


3 ! (m+1 ) !  ,  m+s+2 
(s+m+2) ! 


Inserting  this  weight  factor  in  (10.4)  wo  find  the  approximate  error 

x  +h 

Ri3(x)  ;.'R*g(x)  =  (x 0+h-i)S2j(f  )d^  = 

xo 


(10.5) 
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-  *8+1  ■ 

.s+1  (m+1 ) !  -  Ls+1  i'^8+1 

"h  Tl-Ix,y(*)*-S  y‘-!x.?(x) 


.8+1  (a+1 )  I  (n+1  ) !  f-  !•  i  v  *1 
h  (s+m+2)!  <lYi,s+lJ  -/ _v”| Yi ,s+1 J 


x,y(x) 

(of.  (10.2),  (II.2.3))  or,  ov/ing  to  (II. 2. 5) 
(10.6)  Ri8i^h8+1  i  -IF 


s+1-!xJ(x)t  ‘ 


III.11.  Numerical  Examples 


With  several  simple  examples  having  known  solutions  we  studied  the 
efficiency  of  formula  (2.1 ' )  with  (4.1 ) ,(4. 1 1 ) , (8.2)  and  of  the  re¬ 
mainder  (10.6)  by  means  of  the  subroutines  represented  in  Knapp-Wanner 
/ 29/  or  Y/anner  /51/ •  In  particular,  we  examined  the  question  whether 
increasing  s  and  simultaneously  decreasing  m,  so  that  the  total  order 
m+s+1  remains  constant,  has  a  favorable  effect  on  the  result  or  not. 

In  cloven  out  of  twelve  arbitrarily  chosen  examples,  the  result  was 
positive,  whereas  only  in  one  p.  higher  number  of  Taylor  terms  turned 
out  to  be  more  expedient.  Here  arc  the  results  of  the  example 

l)  y’  =  1-e”y(sin  x  -  cos  x)  ,  y(0)  »  0 

with  the  solution  y(x)  -  log(sin  x  +  ex)  .  The  data  given  are  the 
size  h  of  the  single  step  that  was  calculated,  the  orders  m  and  a  of 
the  formulas  (4»l)»(2.1*),  the  actual  errors  of  the  Taylor  series 
y(x)  with  m  terms,  the  errors  of  the  Lie-series  solution  ( 2 . 1  * ) »  and 
the  estimate  of  the  error  given  by  the  program  according  to  (10.6)« 
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h 

m 

0 

!_ 

Error  in  y 

Error 

in  y 

Estimatos  of  error 

0 , 1 25 

18 

0 

7,?  . 

10-^ 

1  • 

3,1  . 

10-” 

3,1  .  10“17 

13 

5 

2,3  . 

10-11 

1,45. 

O 

CM 

1 

O 

1,0  .  10“20 

j 

8 

:  10 

1 

8,2  . 

10"fl 

4,9  . 

10-21 

1,5  .  10-21 

0 

«o 

*nJ1 

0 

18 

0 

3,25. 

10*9 

2,0  . 

10*11 

1,9  .  10”11 

13 

3,19. 

10’7 

6,6  . 

10-1* 

3,3  .  10’15 

8 

10 

3,62. 

10“3 

1,8  . 

10-'5 

0,2  .  10”15 

— 

; . *• 

-  -• 

- . 

1 

0,500 

16 

0 

1,3  • 

10-5 

6,3  . 

10‘6 

6,2  .  10’6 

13 

:  5  ; 

4,1  • 

10-H 

1 ,1  . 

10’9 

4,7  .  10*10 

8 

10 

1,5  • 

10*  ? 

v»5  • 

10-10 

-10 

0,08.  10 

. i 


Taylor's  scries,  which  converge?  only  for  h  *0,5885...  cannot  bo  user’ 
for  h=0,5  (also  0,25).  Y-.it,  the  Lie- series  correction  gives  good  re¬ 
sells.  Estimation  of  the  error  ie  satisfactory,  especially  in  the  car 
of  t  reasonable  step  size  end  fo?  the  (ujual)  smaller  values  of  s  end 
g. enter  values  of  m  (cf.  Sec.  Ill. 9). 


? '  Oonparison  with  KunKe-Kutta-Fehlbers; 

Differential  equations  of  restricted  three  body  problem: 


1 


—  » T 

*  J2 


y*  =  y1  +  2y.'4  -  * 


(w’=l-y) 

’M  * 

((y^)2*;'")3'2  ’  ( > )2+y|) 37? 


y1+ti 


y1-w’ 


5  =  *4 


=  y  “  ^$2  m§  ^ f 


'j- 


V-, 


((y-j+v)2^?)*72""  7(y ’ 


3-  m 

V<i’*h  this  equation,  in  Durh.ua  /10/  a  comparison  of  different 
methods  was  carried  out  and  there  the  method  of  Fohlbcrg 
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turned  out  to  be  the  best. 

Using  their  initial  values  for  three  different  Arenstorf  orbits 
we  reran  these  examples  with  cur  method,  taking  m=13,  s=3. 

For  comparison  the  results  arc  presented  in  the  following  table 


Table 


The  errors  Ay^  represent  the  amount  by  which  the  initial  con¬ 
ditions  failed  to  be  duplicated  at  the  end  T  of  the  orbit 
(yj-axis  crossing)  for  the  y^  coordinates  respectively. 
Fehlberg's  results  are  taken  from  the  above  mentioned  report. 
Figures  of  the  orbits  and  the  exact  initial  datas  can  also  be 
found  in  /51/„  p. 106-110. 


orbit  number ‘ 


of  1 

positions 

velocities  1 

method 

1 

1 

;  steps 

i 

1  4yJ  ! 

Ay2 

- 

4y„ 

.  1 

1 

“269 

0.3 

|not  given ! 

I0.07 

IT 

1 

Fehlberg 

233 

0.005 

|0.010 

0.004 

0.005 

Lie-series 

2 

395 

0.05  ] 

not  given 

0.1*) 

1 

Fehlberg 

219 

0.009 

0.028 

14.635 

0.381 

l.  . 

'Lic-seno3 

3 

284  | 

0.1 

not  given  j 

i  .783 

i 

2 

i 

|  Fehlberg 

| 

214  j 

0.005 

0.011  j 

0 

• 

o\ 

VO 

Lie-series 

■  i 

in  units  of  10 


*)  These  values  are  probably  not  correct,  since  the  connection 


matrix  for  orbit  2  after  one  period  is 


H(T)  = 


I  3.10*105 
1.72*106 
1.05*10^ 
.  76*10^ 


6.03*10° 

3-l4*103 

1.91O101 

9.23*102 


-9.37*102 

-4.88*105 

-2.97*103 

-1.44«105 


-1.91*10a 
-1.06*104 
-6.50*101 
-2.93*103  / 


This  has  been  calculated  with  using  6  terms  in  each  step.  ThU3 , 
the  derivatives  of  j^CT)  with  respect  to  the  initial  values 
(second  row)  arc  dominant.  During  one  period,  six  digits  arc  1<. 
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3)  Example  for  a  boundary  value  problem: 


j\  =  y2 
y2  =  ya 


(0>  =  y.  ri)  =  0,  y2(0)  =  ? 


:>ineo  y?o  is  not  known,  we  r.ucfis  y2o  and  calculate  the  corres¬ 
ponding  trajectories  y* (x) ,  y2(x).  y*(D  ^  0  we  correct  y2o 
with  Wcwton’s  method 

y'd) 


io(l) 


it' 


24 


Vue  convergence  was  as  follows  (y=10  u  ,  total  time  5  seconds) 

■!?.o 

-0.41 

-0.46358 

-O.46363259167 

'0.4636325917242622617311143 
-0. 46363?59tT2426226l7313495. 

ibiSj  of  course,  is  a  simple  example  only.  Other  examples  for 
bouadury  value  problems  are  carried  out  in  1221  p. 73-94. 
all  computations  wore  carried  out  in  double  precision  (26D)  on 
the  CDG  3600  at  the  Mathematics  Research  Center,  Iladison 
'/is  cons  in . 

■■•'urt her  Examples  can  be  found  in  the  Chapter  on  step  size 
control  ( Ch .  VI.). 
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Chapter  IV 


arCbnor’s  Integral  Equation  and  Convergence 

Proofs 


by  0.  Wanner  and  II.  Reitberger 


In  this  chapter  we  give  a  now  proof  of  the  integral  equation 
of  W.  GrGbncr.  It  is  a  generalization  of  the  well-known 
"variation  of  constants  formula"  to  nonlinear  cases.  It 
makes  possible  an  easy  approach  to  the  formulas  of  the 
preceding  chapter  and  to  a  number  of  further  methods. 

It  also  loads  to  many  iteration  methods,  for  some  of  which 
we  give  convergence  proofs. 

Our  thanks  go  to  Prof.V/.  Crobnor,  K.Ii.  Kast lunger  and  K.  Eglo 
for  their  helpful  discussions.  Wo  further  wish  to  a  c  know  led  p- 
the  suggestions  of  Prof .W.  iiahn,  Grao. 


IV.  1. 
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TV.  1  _The  Integral  Equation  of  Groebnor 


tn  this  equation  appear  derivatives  of  tho  solutions  with  respect  to 
the  initial  values  y^.  Therefore  in  this  chapter  the  following  changed 
notation  is  preforablo: 

- -z.  f(X,x,y)  reap.  Y^A,r,y) 

the  solutions  of  the  differential  equations  (i.  1 .  tresp. (II I>1,i).  h*cn; 

d.1)  f(x,Y(X,x,y)>  f (X,?(X,x,y) ) 

with  tho  initial  values  x,yj  thus  with 


(1.?)  Y(x,x,y)  -  y  Y(x,x,y)  y 

This  moans,  that  tho  dopcndance  of  tho  solutions  on  the  initial  value: 

y  are  now  kept  in  mind.  Special  j  eat  ion  of  these  to  the  prescribed 
initial  values  x..i  y  leads  to  the  functions  of  the  preceding  chapters 

(1.3)  Y(x,xo,yo)  =  y(x),  Y(x,xo>yo)  =  y(x) 

The  connection  between  the  wanted  solution  Y(X,x,y)  and  the  approxinsv; 
solution  Y(X,x,y),  which  is  assumed  to  bo  known,  is  given  by  the  fell: 
wing  thoorem: 

nro  continuous i  then  it  holds  that 


Theorem:  If  f,  f  and  ? 


A 

f  *' 


0.4)  Y(X,x,y)  *  Y(X,x,y)  +  J  f D,  V(X,s,y)'i  ,  di 

x  2  •  •  i ,f(f,*,y>  5 


where 


V  It  • 


Py  'of :  From  c:  C  it  follows  that  l'(X,x,y)  <-C^  (  cf.  /  6/,  p,  2 

j  Y(r  »x,y)  Y(Xl: 


"!c  nov;  differentiate  the  identity 
Y(X,x,y)  =  l(X,w-fY(^,x,y)) 


• » ,y . 


v/i 
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IV. 1. 


with  respect  to  £.  and  after  that  put  £  =x  : 
o  .  *  t(x,y)  dY(a*X'yl 


1  -  A. 

This  is  possible  because  of  Y €  C  .  Finally  we  insert  £  for  x  and 

•> 

for  y: 


A  similar  differentiation  of  Y(X,^,Y('; ,x,y))  with  respect  to  {E  yield.’ 
(using  chain  rule  again) 


Y(Xff,Y(*.,x,y)) 

v 


Finally  we  subtraot  this  from 


»U*.x.y )  +  J(  v) 

a,  tf(x’y)  ay  Jt,?(|l3[>y) 

(1.5)  and  integrate  from  x  to  X: 


*  [-Y(X,|,Y(|,x,y))] 

-  -Y(X,x,y)  +  Y(Xfx,y) 


-  Y(X,X,?(X,x,y))  +  Y(X,x,Y(x,x,y)) 
(cf.  (1.2)). 


Thus,  (1.4)  is  proved.  The  different  arguments  £  and  x  in  v'Xf£,y) 
do  not  mind,  since  they  are  equalised  by  the  substitution  rule  x  ->V> 
y  ->Y(fc,x,y). 

J  Done. 

This  integral  equation  was  found  first  in  i960  by  Groobner  for  ana¬ 
lytic  equations.  It  was  rediscovered  in  similar  form  (cf.(5.l))  in 
1961  by  Alekseev  / 5*1/  •  The  above  given  proof  is  similar  to  that  of 
Alekseev. 


IV. 2  A  Generalization 


The  above  integral  equation  can  be  generalised  in  the  following  way: 
Theorem:  If  f.  are  continuous,  F(x,y)  is  oontinuously  diff erenti-ib  1 

it  holds  that 


F(X,Y(X,xfy))  =  P(X,Y(X,x,y))  + 

X 

(2,)  ♦  J[b2  F(X,Y(X,-,y))j  ^  ^  aJ  • 

J  , 

I 

Clearly  (2.1)  ooinoides  with  (l.  J  if  F(x,y)  a  y.  For  analytic  equations 
this  formula  has  first  been  recognised  by  K.  Egle  (  of.  QSR  Nr.l). 

Proof:  First  differentiate  F(X,Y(X,>:,Y(f ,x,y)))  with  respect  to  £ . 


|^P(X,Y(X,|,Y(^,x,y)))  = 

/9F(X,Y(X,£,Y(£,x,y))).  faY(X,_x„y)  +  f (x,y) 
ay  *'  1.  rx 


Pext  we  multiply  (1.5)  by  a— •  F(X,  Y(X,  £,Y(|,x,y)) )  and  subtract  the  tr.c 
formulas  : 


-  F(X,Y(X,|,Y(^,x,y)))  = 

-  j(f(x,y)-  f(x,y))  ~U,Y(X,x,y)>  aY(^x-«y)  ] 

y  y  Jf,Y(f,x,y) 

=  l(f(x,y)-  f(x,y))  P(X,Y(X,x,y))l  A 

"  dy  -?,Y(vt,x,y)  . 


>  V‘ 


Integration  from  x  to  X  now  yields 
3inco  again  x 


[-  F(X,Y(X,j,Y(^,x,y)))]  - 

+  F(X,Y(X,x,Y(x,x,y) ))  -  - 


the  wanted  integral  equation  (2.1), 

-  F(X,Y(X,X,Y(X,x,y)))  + 
(X,Y(X,x,y))  +  F(X,Y(x,x,y)) 


Done. 


IV. 3  A  Volterra  Integral  Equation 


Interchange  Y  and  Y  in  (1.4): 


This  yields  the  fol. owing  Volterra  integral  equation  for  the  solution 
Theorem  If  f.  f  and  are  continuous,  it  holds  that 


(3.1) 


Y(X,x,y) 


-  Y(X,x,y) 


In  addition,  if  F(x,y)  is  a  continuously  differentiable  function,  wo 
have 

(3.2)  P(x,Y(X,xfy))  -  F(X,Y(X,x,y))  + 


+ 


[D2  F(X,Y(X,j,y))] 


These  formulas  differ  from  (1.4)  and  (2.1 )  only  by  the  exchange  of  f 
* 

and  Y  under  the  integral  sign. 


IV. 4  The  Variation  of  Constants-Formula  as  Special  Case. 


^et  f(x,y)  be  linear  in  ys 


?■ 


A  .  A 


and  let  F(x)  *  (F^^(x))  be  tho  fundamental  system  of  solutions  with 
F(x)  =  I  (identity  matrix).  Then  Y(X,Jfy)  =  F(x)F“1(t)y  ,and 

”  \. 

0 


§£  Y(X,|,y)  -  F(X)  F-1  ( f-) 


Thus  for  the  solution  of 


TV.!?. 
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(3*1)  rends  ns  follows: 

f 4 •  1 )  Y(X ,x . y)  «  F(x)y  +  f  t5’1  (  M  r(fc,Y(f  ,x,y))df  . 

This,  however,  i3  nothin g  else  than  the  Variation  of  Conotants-Pornu] 
for  inhomogeneous  linear  <lif  feiencir.i  ays tons . 

This  shows,  that  (3*1 )  acts  the  sane  part  for  nonlinear  equations, 

than  (4.1)  does  for  lima-  equations.  I.e.,  it  has  applications  to 

asymptotic  theory  of  differential  -equations  (  e.g.  <1  as ov<  /5?/»P*6f»' 

to  s  c  ibility  theory  (P.  Br.auor  I'yj'l ,  nl  ekcoev  / 54/)  ->r  to  the  treat- 

mont  of  stiff  differential  aquations.  This  wo  hope  to  discuas  in  a 

# 

later  report. 


IV, 5  Proof  of  the  Formulas  of  Section  III.? 


The  series  of  the  preceding  chapter  are  now  simply  obtained  by  in¬ 
serting  the  Taylor  series  of  the  solution 


/  r  *  \  V 


v r-.<, ••  j  j  {V”Vi  ci i ) 

“■°  ,  |  rt 

into  the  right  hand  side  of  (1.4)? 


/ 

f?  rt 


(5.?)  Y(X,x,y)  *  ?(X,x,y)  +  ?'  J  fl)0Day]  a  d« 

«  -  C,Y(h>x,y)  * 


+  Vx»x*y> 


n  ■  r 


A 


with  „  ,v  _.x  !  ■;  (X- »',}*'  r„s+1  1  *T 

«  U.,x,y;  -  |  j  !  L..,  -~yu-  U>  y,  .  ,v  d».i  . 

*  r  ■  -  u 


or  by  interchanging  the  order  of  integration; 

X  ’1 


», (x,*,y)  -  )  j  I'd,  iV+V 

J  v-  ».✓  c. 


\ 

il  i  • 


x  x 


!  L  *,  gj 


»Y(vj,^,y)  k,Y(£,x,y)' 


The  inner  integral  of  this  c^n  no  -  he  evaluated  with  the  help  of  «■'. 
vnoralizod  integral  .  qua  ti or  (?.')  to  give 
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Theorem:  If  f ,  ?,  ~  are  continuous  and  P8+^y  is  continuously  diffe¬ 
rent  iable.  formula  (5.2)  is  valid  with  the  remainder 


(5.3)  Ra(X,x,y) 


I  ,  -b 

x  *■  -.],Y(Tj,x,y) 


s+1y]  a  ]  dt| 


Proof :  Under  the  given  conditions  (5*1)  and  (1.4)  are  valid.  (2.1)  is 
used  Tilth  the  funotion 

b+1 

F(x,y)  ■  1  g  |  D  y  and  with  X  roplaoed  by  VJ.  The  stated  c> 

ditions  allow  this  application.  Done 

(5.2)  and  (5*3)  are  nothing  elso  than  the  formulas  of  III.3»  if  trans¬ 
cribed  to  the  original  notation. 

Remark :  Knapp ,/26/,  has  proved  (5*3)  under  the  weaker  condition  D  J*  C 
Following  we  prove  the  formula  of -Groobnor  in  III • 1 : 


IV. 6  Convergence  for  s  a> 


Theorem:  If  the  functions  fi(x,y)>  i.e. .  the  operator  D,  are  analytic 
in  some  domain,  then  for  a  sufficiently  Bmall  h-X-x  (5*3)  converges  to 
zero,  i.e.  we  have  from  (5*2) 

o-*  X 

(6.1)  Y(X,x,y)  -  ?(X,x,y)  +  2.  J  -^^fD  Day]  ^  d£, 

«-o  x  a!  L  Y(^  ,x,y)  ' 

the  formula  of  Groebner  stated  in  III.1. 


Proof :  Since  D  is  analytic,  it  can  bo  majorized  by  the  operator  whi 
is  in  one  variable  only 


I)<  & 


N  d_ 

0-  f) ,,z 


(of.  o.g,  Groebner  / 22 /,  p.  30jand  Groebner-Watzlawek  /22 /,  p.  22  5). 
Thus 


IV. 6. 
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Remarks ; 

1)  A  second  proof  of  (6.2)  is  possible  by  inserting  the  infinite 
power  series  into  the  integral  and  assuring  uniform  convergence 
.rhich  allows  interchange  of  summation,  integration  and  differen  ■ 
ifttion . 

,  )  Jtill  another  proof  (the^  historically  first  one)  was  given  bp 
tirdbner  by  rearranging  the  power  series  for  the  solution  Y  in  a 
special  way  (of.  e.g.  Grfibrer  /22/,p.35,  Knapp-Wanner  /29/,p.29). 


V 


5? 


IV. 7. 


jlV.7.  A  General  Process 


Formula  (5.2)  has  resulted  from  inserting  a  Taylor  series  solu¬ 
tion  into  the  integral  of  Grdbner's  integral  equation  (1.4). 
This  leads  to  the  idea  to  insert  any  approximate  solution 
Y(X,^y),  say,  of  order  s.  We  thus  obtain  a  new  approximate 
solution  Y(X,x,y)  given  by  the  formula 


X.  * 


(7.1)  Y(X,x,y)  =  Y(X,x,y)  +  /  jp2Y(X,e,y)] 


K >Y(C,x,y) 


dt 


A 

Theorem.  If  Y  is  of  order  m,  Y  is  of  order  s,  then,  under 
appropriate  differentiability  conditions.  Y  is  of  order  m+s+1. 


Proof : 
to 


rst 

Because  of  the  order-condition,  the  error  of  Y  is  equal 


error  of  Y(X,^ ,y) 


s+1 


(X-^) 
TiTiTT 


F(X,^y). 


The  error  of  Y  is  now  obtained  by  subtraction  of  (7.1)  from 
(1.4) 

X  ^ 

error  of  Y(X,x;y)  =  /  [b3(  error  of  Y(X,e,y))J  ^  d£. 

x  C,Y(e,x,y) 


Again ,  as  in  section  III. 6,  contains  the  factor  U-x) 
and  wo  thu3  have 


m 


error  of  Y(X,x,y)  =  /’  0(e,x,y)  [|^F(X,C ,y)]  . 

Since  (X-Oa+^(C“x)m  does  not  change  sign  in  the  integration 
interval,  the  moan  value  theorem  can  be  applied  and  yields 

X/V  _\S+l/_  vm 


error 


?(X,x,y)  =  /  G(e,x;y)[|-P(X,0..y)j  f 

x  8  f  Y 


(X-x) 


m+s+2 


lm+s+2) ! 

thus ;  Y  is  of  order  m+s+1. 


(x<0<X) 


Done. 


IV. 8. 
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Remark .  The  Theorem  in  section  III. 5  is  related  to  this. 

Hence ,  to  each  pair  of  methods  with  orders  m  and  s  and  with 
solutions  Y  and  Y  resp.,  formula  (7.1)  leads  to  a  new  method 
with  order  m+s+l  and  solution  7.  The  Lie  series  of  Chapter  III 
is  obtained  by  inserting  m  and  s  terms  of  the  power  aeries 
expansion. 


IV. 8.  Iterated  Integral  Equations 


Still  further  integral  equations  are  derived  from  (l.*l)  and 
(2.1)  by  iteration: 

-  X  _  . 

(0.1)  Y(X,50,y)  =  Y(X,C0,y)  ♦  /  LDjYCX.q.jrV.L  + 


X  X 


*1  /  [°2[V:(X^2*y)^  -L  as2dsl  - 


«o  «1 


Z2  Z1 


and 


(3.2)  F(X,Y(X,s0,y))  =  F(X,Y(X,c0,y))  +  /  [t^FtX.YU.e^y) )]  a 


X  X 

♦  /  /  [D2[D.P(X,Y(X,c2,y))^  ds^j 

'o  51  z2  Z1 

where  z^Zj***  denote  the  following  substitution  rule 
( 3  •  j )  zk  — 


Still  more  general  equations  are  derived  after  repeated 
iterations. 

The  following  theorem  is  obtained  by  the  repeated  application 
of  the  preceding  theorem  (section  IV. 7): 
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Theorem:  If  Y  is  of  order  m,  Y  is  of  order  s,  then,  under 
expropriate  differentiability  conditions,  7  which  is  defined  l>y 

■•r 

(8.4)  Y(X,e  *y)  =  Y(X,c^y)  +  /A[D2?(X,c  ,y)]^  dc  ♦ 

-  5o  Z1 

♦  /“  JX[D2rP2Y(X,5  y)]^  dc2dsl 
*o  «l  z2  Z1 

1;  of  order  2m+s+2.  Done. 

Voirgrks :  1)  The  insertion  of  a  power  series  for  Y  leads  to  the 
series  with  multiple  integrals  as  given  in  Wanner  /51/,P*73“74 

2)  The  order  of  the  analoguous  formula  with  an  r-fold 
integral  is  rm»s+r . 

IV .J_-  Iteration  methods  and  convergence  proofs 


The  formulas  (7.1), (8. 4)  etc.  can  be  iterated  in  many  ways. 

I  acre  are  possible  iterations  with  respect  to  Y,  or  with  respect 
to  Y,  cr  both.  A  number  of  methods  appear  as  special  cases .  such 
ue  the  iteration  methods  of  Picard  or  that  of  GrSbner-Knapp ,  the 
method  of  Poincare  and  so  on.  Por  a  few  of  these  iteration 
methods  convergence  proofs  and  error  estimates  now  are  given, 
for  others  we  have  not  yet  found  them. 


IV . 10.  The  Iteration  Method  of  Grobner-Knapp 


This  method  appears,  when  (7.1)  is  iterated  with  respect  to 
Y  while  for  Y  the  first  s  terms  of  the  power  series  solution 
are  insex’ted.  Thus  we  have  the  iteration  process  (cf.(5«2)) 


(tO.l) 

Y:r^1)(X,x,y)  =  Y(r) (X,x,y) 


s 

l 

1=0 


L(x-c)‘ 


[P 


(r)^a 


o“y]  (r)^ 


IV. 10. 
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i  - 


Y(r)'  .  J(r)(x//(r))#  0(r)  s  |_  4  f<r)<x4y)  |-  , 


starting  with  Y^°^  =  Y. 


The  convergence  follows  from  rouation  (III. 5. 3)  under  the 

condition  that  D‘  y  satisfies  a  Lipschitacondition.  If  the 

Lipschits-condition  is  assured  only  in  some  compact  domain  B 

,  * 

(as  usual  v/ith  nonlinear  equations) ,  further  considerations 

tr\ 

are  necessary  to  assure  that  the  iterated  functions  Yv  '  do 

not:  leave  B: 

< 

Theorem:  Assume  that  the  Lipachitzcondition  (111.5*1)  for  the 

+  1  ~  -  i— rrr 

functions  D  y  is  satisfied  in  the  domain 
(10.2)  B  =  {(5,r.)  |x«£sx+a,  ly^-n^lsb) 

A 

and  that  tne  approximate  solution  Y  satisfies 


CIO. 3).  i^U.x.y)  Y.(c,xsy)|  s  H 


m+1 


nnr 


n  n  .  o 
“  (ri+1 ) !  *  72 


for  X&C4X+&.  Further  h-'C-x  she. aid  satisfy 

,,rti+ 1 

(10.  '4) 

(10.5) 

(10.6) 


!?;  (e/-.y)  '  y- !  *:  f 

A  %  (L 

K3nhs+1 

TnH^(m+T)Tr.Tit+s7^  '  1 

(r). 


x<f,.5X+h 


Then  the  iterated  solutions  Y  y(t,x,y)  of  (10.1)  do  not  lcav 
B  for  x^t<>:+h  end  the  iter  tlons  are  arbitrary  often  possible 
and_converge  to  the  solution  with  the  error  estimation 


(10.7)  n^U.x.y)  -  Y-r)U,x  „)j  5.  M!c-*rT1  fpil .TSTf 5F5T  • 


m+l  tVI«-xi,+1>r 


free: 


•  C  y*  } 

.)  First  we  snow  that  the  functions  Yv  ’  are  again 


(c+1) -times  differentiable ;  since  only  for  f the 


formulas  (5. 2), (5. 3)  are  valid.  This  wo  3imply  show  by  (s+1)- 
times  differentiating  Rg  in  (5*3) s 
Using  the  well  known  formula 


,s+l 


dX 


s+T 


X(X-e)3 
s ! 


B(€)<U  =  B(X) 


we  obtain 


,3+1 


^■571  s  fCa+ly]XjY 


[Ds+1y] 


x.y 


which  is  continuous.  Now  the  assertion  follows  from  the  fact 
that  Y(X,x,y)c  CB+1(cf.  section  II. 1). 

2)  Next  we  confirm  that  the  iterated  functions  do  not  leave  B: 
First  by  (III. 5. 3) ,(10. 4), (10. 6) 


m  M|c-x! 

(10.8) 


m+1  K  n|5-x|s+1 


(m+l)l  (m+2) . . (m+s+2) 


s  |  i 


4 


si 


again  by  (III. 5. 3) 
(10.9)  |Yi-Y^2)|s 


II  |  C-x  | m+1  Kt.n|5-xjs+1 


Ksn|(-x 


s+1 


(m+1)!  (m+2) .. (m+s+2)  (m+s+3) . . (m+23+3) 

' - ' -  ' - • - ' 

J  <i  <i 


and  so  on. 

Thus  by  (10,5)  and  the  triangle  inequality 


*yryi 


(r) 


I  *  iyrYil  + 


4 


y(r) 

Yi 


*  b. 


(10.7)  follows  by  induction  from  (10.8) , (10.9) ,. .  by  (III. 5. 3). 

Done. 

Remarks .  1)  The  number  m  in  (1C. 3), (10.4)  need  not  be  the 
greatest  possible.  Perhaps  sometimes  (10. 4)  may  be  less  re¬ 
strictive  for  smaller  m;  (10.6)  however  is  not. 


Mo” 


IV. 12. 
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2)  The  theorem  is  essentially  due  to  K. Knapp,  the 
•.  r-oc f  is  new. 

3)  We  already  have  a  convergence  proof  for  an 
arbitrary  Y,  not  only  for  finite  sections  of  a  power  series, 
and  we  are  at  the  present  working  on  its  simplification. 

» 

IV. 11.  Picard’s  Method  as  Special  Case 

* 


The  method  of  Picard  comes  out  from  (10.1)  by  putting  s=0: 


Y'*T*'(X,x,y)  =  Y(r) (X,x,y)  +  /  [Dir;y]  (r%dC 

x  f,Y^  7 


V 


,(r), 


J*  f(c,Y(r))dC-  /Xf(r)(C,Y(r))de 

>:  x 


X 


Y(r)(X)+y 


=  y  +  /  f(5,Y(r'(C,x,y))d€, 


ihicli  is  the  well-known  iteration  of  Picard., 


tV.12.  Poincare's  iicthod  of  Parameter  Expansion 


I 3  (71)  iterated  with  respect  to  Y  while  Y  is  kept  fix,  we  gat 
Poincare's  method  of  parameter  expansion;  i.e.,  we  obtain  the 
solution  expanded  in  powers  of  a  small  parameter. 

To  show  this,  assume  that  the  operator  is  multiplied  with 


a  small  parameter,  say  e 


D2 


«  sD. 


# 
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X 

(12.1)  Y(r+1)(X,x,y)  =  Y(X,x,y)  *+  /  [D*Y(r)(X,5,y)]  <U 

x  €  ,  x 

starting,  say,  with  Y^°'sY.  Inserting  for  Y^  again  and  again, 
it  can  be  seen,  that  Y  is  expanded  in  powers  of  e(in  simplified 
notation) : 

(12.2)  Y(r+1)  =  Y  ♦  c/D*Y  +  ^/Jd'd'y  +  ...  ♦  er+1J .  ./D* .  .D*Y  . 


r+1  r+1 


Subtracting  this  from  (G.l)  (more  precisely:  from  the  (r+1) -time 
iterated  equation  similar  *o  (3.1) 

(12.3)  Y  =  Y  +e/DjY  +  ...  +  Er+1/.  ./D^.D*?  +  er+2J. . JdJ. .D*Y  ) 


r+1  r+1 


r+2  r+2 


we  obtain  for  the  error 
(12.il) 


Y  -  Y(r+1)  =  er+2/ . ,  /ol^DgY  . 

r+2*  r+2- 


Wo  again  consider  the  example 

y*  *  /x  ♦  Sy,  y(0)=0 


which  is  known  as  an  equation  which  poses  difficulties  to  a 
numerical  integration  (cf.  Rosser/ll2/,  Cooper/8/).  Again  let 

f  =  /5c,  D2  =  D*  .  /y  ,  t  .  1. 

Then  the  above  expansion  becomes 


y(x)  =  | 


X 


3/2 


(2  4 

hr 


,7M 


4 


f! 


l 


,9/** 


which,  for  small  x,  gives  an  exactitude  of  hundreds  of  Runge- 
Kutta  steps. 

A  convergence  proof  of  series  (12.2)  is  given  in  section  IV.  1^. 


IV. 13. 
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The  actual  determination  of  the  aeries  (12.2)  is  mostly  easier 
by  the  usual  expansion  as  e.g.  it  is  described  in  Knapp- 
banner  /28/,  section  IV. 2.. 


IV. 13.  Power  series  as  special  case 


Put 


'i(x*y)  sfp 


then 


Y(X,x*y)  =  ya 


and  the  process  (12.1)  reads  as  follows 

(13.1)  Y(r+1)(X,x,y)  =  y  +  /Xf(C, y)  ~Y(r)(X,e,y)d? 


a' method,  which  is  connected  with  a  power  3erirs  expansion 
in  the  following  way : 

Proposition:  In  thw  case  cf  autonomous  equations.  (13.1) 

coincides  with  the  method  cf  power  series: 

This  is  seen  by  induction:  For  autonomous  equations  f(x*y)  does 
not  depend  on  x  and  we  have  for  (13*1) 


y(r+!)  s  y  +  f(y)  jX  Y(r)(X,F,>y)dc. 

*' — v— 1 *■*'  X 

D 


Thus  if 


then 


=  V  pay 

.So  al 

y(r+1^  .  y  ♦  I  D  /  -  d5  Day 

01  =  0  X _ ° _ _ : 

~’7x-x)a+1 

hottt 


and  with  0=a+l 


IV. 14.  Convergence  proof  of  Poincare's  method 


The  following  convergence  proof  of  the  iteration  (12.1)  a 
starting  with  an  arbitrary  analytic  function 

X(o)(X,x,y)  =  l  i~Ti“5Jy 
v=o 

has  been  worked  out  together  with  K.  Kulmert . 

First  get  a  majorization  operator 


A  = 


N  d_ 
(i-z/p )  dz 


for  Ds  3  as  well  as  cD2  =  Dg. 

Then 

(14.1)  Y(X,x,y)  -  Y(r)(X,x,y>  -<  2  J  H ^-X-gP-Aaz . 


where  the  symbol  -<  denotes  rasterization. 

(14.1)  is  proved  by  induction: 
r=0 : 

Y(X,::,y)  -  Y(o)(X,x,y)  =  J  (Daz  -  U® z) 

/M  • 


a  =  0 


■< s  L  ^  • 


y=o 


=  1 


r+r+1 : 


Y(X,xJy)  -  Y(r+1)(X,x,y)  =  ej  [D5(Y-Y(r)JA  dc 
X 


\ 


A  ^  /  V  v  ^  01 

/  [42  I  (r)^~ 
x  a=r 


E,  I 


y-c 


Y  J 


1 

d$  . 


Application  of  the  comutation  theorem  (cf.  e.g.  /22/,p,17) 
and  rearrangement  of  the  double  3um  gives 


W  t 


4 
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<2  i  j  r»i  ;x iyai d5.4^“nz 

.i  t.!  y! 


a=r  y=o  oi 


(X-x) 


( a+Y  +  1 ) ! 


and  with  3=u+1+y 


i-Y!r+1)  -<  2  I  T/")1— Sr  A6z 

gsr+1  &=rvrv  y 


(r+l) 


Done . 


Next  wo  transform  the  initial  values  y  to  the  origine  and  sum 
up  (14.1)  for  rsO,l,..: 


“  “  /  ,,  x  u 

I  I  (“ 


r=o  a=r 


S  I  II) 


a=o  (r=c 


which  converges  for 


|x-x|<^j  . 


Thus  it  is  necessarily 

lin  (Y-Y(r))  r  o. 
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Chapter  V 


Runge-Kutta  Processes  with  Multiple  Nodes 


by  K.H.  Kast lunger 


Abstraot ; 

®he  use  of  the  differential  operator  D  makes  it  possible  to  extend 
the  methou  of  Rungo-Kutta  in  suoh  a  way,  that  the  power-series 
method,  the  classical  Runge-Kutta-method  as  well  as  the  proces¬ 
ses  of  Fohlberg  are  contained  as  speoial  oases.  The  generali¬ 
zation  is  in  such  a  way,  that  not  only  the  functions  f^Xjy) 

are  evaluated  at  some  intermediate  points, but  in  addition  also 

2  m 

the  functions  ,  D  f^  ,  ,  35  f^  .These  new  methods  are 

advantageous  especially  when  combined  with  the  conoept  of  re¬ 
cursive  generation  of  the  values  of  35nf^  »  35n+Vi  »  as  desribed 
in  Chapter  II  . 

We  first  develop,  the  general  form  of  the  conditions  for  the 
coefficients  of  those  processes  thereby  extending  the  results 
of  J.C.  Butohor  /  1/,  These  equations  become  still  more  compli¬ 
cated  than  those  for  classical  processes. 

Next  it  is  shown  that  to  each  quadrature  formula  with  multiple 
nodes  there  exists  an  analoguous  Runge-Kutta  prooess  with  the 
same  number  of  nodes  and  with  the  same  order. This  again  extends 
results  of  Butoher  /  2 / . 

Fohlberg1 s  method  is  shown  to  bo  nothing  else  than  a  generalized 
Rungo-Kutta  process  with  one  m-fold  first  node  and  a  few  additio¬ 
nal  single  ones. 

Finally  we  give  examples  of  expllcite  prooess  and  numerical  example  . 


V.l. 
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V.l.  General  Theory 


Notation 


Let  the  following  autonomous  system  of  ordinary  differential  equation? 
he  given 

(1.1)  yj-  fi(y1,...,yn)  or  for  short  y'»  f(y) 
with^tho  initial  conditions 

(1.2)  ?!(*„)-  Vi0  Ol  y(x0).  yQ  . 

Here  we  treat  autonomous  systems,  since  then  the  following  theory 
becomes  more  simple  .This*,  of  course,  is  no  reduction  of  generality, 
what  can  bo  seen  by  introducing  x«y  as  new  variable  with  y'  -  1. 

Are  the  functions  f^  analytic  in  a  neighbourhood  of  yiQ,  then  the  solutir  i 
of  (1.1,2)  is  given  by  the  following  power  series 


(1.3)  yOO* 


k»o 


o 


7ith  the  differential  operator 
n  a 

T)  "  8 

fa >  *3 


(1.4) 


Again  the  symbol  [, ..]  means  that  after  all  differentiations  the  initial 

values  x  ,y  arc  to  be  inserted, 
o  o 


Elementary  Differential s 


Definition: (Butcher  /  1  /,  p *  1 87 ) 

{1}  :=■  f  is  the  only  one  elementary  differential  of  order  one; 


(1.5)  F«  {?,... Pe  )  1 


n  .s 

i?  p  P, 

- rrr  ‘'V  SJS 


J1 

where  F^  (Fi1 , . . . ,Fin) 

is  an  elementary  differential  of  order  r  and  degree  8  , if  F^  are 
elementary  differentials  of  order  r^  : 


order  r:  r«=  r.+...+r,+1 

•  —  IS 


degree  s:  8=  number  of  F^,  which  constitute  F  , 


v 


more  generally: 


(1.6a)  P-{P/...P  ®)  { F4 .  •  .F,.  •  •  .F  . .  .F  ) 

1  °  vJL^ _ h  kSL^, \Sj 


H1  *0 

Is  a  elementary  differential  of  order  r  and  degree  s, 
whoro 

(1.6b)  r«  ^r^. .  •+nffr(J+1  and  a-  . 

In  this  notation  the  exponent  may  be  zero  also  j  we  also  define 
{F°...F°}  s-  f  . 

Notation  convention: 

For  simplification  we  now  introduoe  the  following  notation: 

if  v1f...,ve  are  vectors  v^-  (v^ , . . . »vin)»  'then  we  denote: 


(1.7) 


d°  n  da 

r-#.tV  -r—  y*  V.  .  .  .#V  .  -r- - 5“  I 

ay  i7T....js-i  3i  83b  3ya,"' 


more  generally:  with 

a 

(1.8a) 

we  write 

(1.8b)  v^...v^  ^7  :* 


iV:;vi;  * 


iToing  this  notation  (1.5)  and  (1.6)  now  becomes 

I*.  4—*  • 


d8f 


0.9)  {V..P.|.  . 


.  lp, ••••»„' )  ■  V— p 


^  dSf 

<3  .  S 

dy 


The  following  theorem  about  elementary  differentials  is  due  to  Butcher  /  1 

Theorem:  D  f  is  a  linear  combination  of  all  elementary  differentials 
F  of  order  r  with  positive  integer  coefficients  a: 

(l.lo)  Dr-1f«  X _ «F  *  ^ 

Ord  F-r 

The  ooeffiolent  a  of  P-  {  F1  .  ..F^  )  is  given  by 


1) 


The  symbol 


denotes, 


that  the  summation  is  over  all  elementary 


differentials  of  order  r 
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(1.11) 


*i 


<r-1>!  J  Liri  ferJ  • 

1=1  1  i  \  i/ 

where  are  the  coof  icients  of  their  orders 

and  r  the  order  of  F. 


Examples : 

Df=  {f > 

D2?=  {{f  }}  +  {f2} 

D5f=  {{{  f  }}}+  3{{  f  )f>  +{{f2}}+{f5)  . 

A  Xi.it  of  all  elementary  differentials  up  tc  order  eight  is  given 
in  /  I  /,  pp.  191-193. 


next  modify  the  above  theorem: 

p  m  1  p  _ 

with  D  1  f=  D  iy=  ]T _ Oj  F. 


Ord  F^r^ 


i  i 


we  have 


1.1.12)  {  (D  y)...(D  3y)>  =  2 


Ord  F1=r1 


-  •  •  *  ^ 


Ord  F  »r 
8  s 


v--\  (P1--V 


Theorem :  It  holds  that 
r-1 


(1.13)  3)  y»  L  2 


_ T'-“T7^(  (Ar). ..(»%)> 

t-1  r<+.„+r^»r-1  Tl  -r**#V 


1  _ 'r-1 ) ! 


1  ■  •-  '•*+' 
1  t 


Proof :  Inserting  (1.12)  into  (1.1 3)  wo  obtain 

^  ^..-Cr-“lJ-!  .  -J-  X _ *  •  •  /I _ “i  *  •  *at  {  Fi  •  •  * *  • 

t-1  r1+...+rt-r-1  J--d  F1»r1  Ord  Ft-rt 

1 

It  is  easy  soon  that  all  elementary  d ifferontial3  of  order  r  appear 
ir.  this  expansion.  If  we  pick  out  one  of  these,  say  {F^1 . .  .F^0)  , 
v,e  see  that  the  coefficient  of  this  is  equal  to 


i*5lli  -ir  Ctfi 


JJ,  PJ J  =  (r-,)! 

This,  however  is  exaotly  the  coefficient  a  • 


nfr  i- 
jJi  “j 


!i  “1 

!  '^y 


! 


1.'. 


Power  Series  for  the  Ru ngo-Kutta  Approximation 


The  olassioal  Ruhge-Kutta*method  uses  the  following  formulas  for  tne 
approximation  y(x) : 

n 

(1.14a)  y(x>«  yQ+  h 

and 

(1.14b)  f(yQ+  h  ^biJgj)  . 

These  formulas  arc  now  generalized  in  the  following  way: 

n  »  number  of  nodes (stages)  of  the  method 
m  *  their  multiplicity 

d.i5)  *<*>. 

(1.16)  5(x). 

Remarks : 

1)  For  m«1  wo  get  (1.14)  again. 

2)  For  simplifioity  we  put 


(1.17) 


^*o  for  k-m+1  >m+2 . . .  *  >  i»1  >2 . .  • .  .n 

bij^“°  k»m+1  ,m+2f .  •  •  ,  i,  j»1 ,2 , . . .  fn 


3)  The  assumption  that  all  nodes  have  the  same  multiplicity  is  not 
restrictive.  Otherwise  we  put 

m-  max  {m^,...,mn}  ,  m^«  multiplicity  of  the  i-th  nodo 

(k} 

o^  '*o  for  k*  m^+1 »m^+2 , . . .  ,  i=1f...»n 

(1.18>  (k) 

b^ j  *o  for  k«m jtl ^m j+2 1 • » *  $  i>jB1t»..>n 


Theorem  (Expansion  Theorem  for  Runge-Kutta  ) : 
a)  posesses  the  following  power  series 

(119)  »oo.  ®  si  Roo 


v.1. 


the  vectors  R>  '  arc  determined  rcoursivly ; 
( 1 .2oa)  H ^  -  [Dky]o  k-1,2,„. 


(l.2ob)  R^'  -  y  / _ 7^77  S  ...Si  H  —  ] 

1'U  T>  x1+,m+*1.|;,1".,,xt'  3’KV  i,HT  dyT 

*  *  1  T 

i  k= 1 , 2 , . . .  , 

o-2»  si,„- 

b)  Tho  approximation  of  the  r.unge-Kuttn-mothod  has  the  following 
expansion : 

(1.22)  J(x)-  y  +  TV 


t.  g—j  (Tt-klT  f-j  i  i ,  H-k 


Proof : 


1)  Since  f  is  analytic,  it  follow.**  that  ^  posesses  a  power-series 
expansion  (1.19)  • 

2>  vj,  TT  (’  ,?)  V  i,7T  %'iW  (-‘,9) 

<**>  n  /  -*cc  v,  ji  /  \ 

j  t)  e  j  bj^y  £7  bW  . 
r  £ V-  £i  13  tfeo** !  J'*1 

y  >f  y  f  *<;>»<’>  <■*«> 

n  £l  £o  **,u'!  l!|:!  W—>  Z'v 

<X,M  H  n  (  \  (  \  2  v,H 

V  +  ~~  (M|V  b  T  R  ^  -  7  j-  A—  g  K  y  +v 

yo  “T|H!  i.j  Rj,K-T  '  n'^TiH!  i»K  0  1 


with 


ce  ,  K 
n 


^•24)  vr  ^  —  SUl< 

fl.n)  5.  .  7  (“')7b(T)B,:T) 

1,11  *T/  jtl  1,1  J>H~r 

d _ /  T-  n  ^  jL. 

Tl^*lfrr  -*.V  a7“ 


3) 


•  ' 


V.  1 


»  h 1 


„  (H.+...+X  =|i) 

11-21  ;—•  »■■-  si ,  -Si  „  —z  -  T 
CT  hTT  V***V  1,h1  1,mt  dyT 

*  T 


(1.25) 


oo 


hH 

;4! 


Z. 


uL 


- p‘  ,,  - — - -  x  '...w  !  i,X  ,  T 

U«T  X.+. .  .+K  «(i  1  T  1  *  Tdy 

x|>  1  T 

•I;  Taylor's  theorem  for  multiple  variables  with  use  of  (1.Y  )  and  (1.8) 


;i.2S)  f(vv>-  fi.  ii  v1^  f 

v*  T1  idy" 


;)  e'lL'*”  (D'ty)(y0+v1) 

[*]..  i,  a,  ^ 


T® 

5»  i 


(1.26) 


llDkyl  Oi*« 


dyT 


[*].♦££,££  2L_ 

u-r  k,+.  . . 


JiL 


- - ;t1 1...H  !  SifK1,,,Si,H 

1T...+K  =(1  1  T  1  T 

x]  »  1  T 


r  T 

s* 

ayT 


Jo 


[iky] 


Ct> 


..  ...S,  ..  Dky 


0+  ^TT  ^  v-T  ”!  x(‘+7.  '."+x  =  MH1 !  ^ MT 1  Si » H1*  *  ‘ Si ’ Xt[dyT 

Hi  1  T 


!  0 


A  comparison  with  (1.19)  gives 
(k)  nk 


(I*20®)  *i,„-  [*]„ 

(i.?ut)  r.W.  T  1,  T 

1,11  vrTl  £f777 


i 


_ u  i  .  ,  S,  ,  ...s,  y 

-n  *i  ••••%•  1>H1 

1  T 


(dy 


o  yw- y0+i>y:«|1)41)+"-+i'"iioim)«im)  (1i17) 

i=ik)4k)  (1^19) 

v&VV0*#  - 

k«1  i- 1  li=o  r 

..  .V-  hk+ti  (k+tx) !  T^_(k)„(k)  (k+iiex) 

U  ^  6rl  ' 

y  +7" ill  hi  y“cOORv<)  „  y  +y  ^  T 
r  £tbi  frr 1  i>x-l!  0  »iK!  H 


.1. 
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with 

(1.23) 


T  a  r~  .  2tl-  y-0(k)R(k) 

Th  ^T^Syr^ri  Ri,x-k  * 


Done. 


Connection  w ith  Elementary  Differentials 


Theorem:  is  a  linear  combination  of  all  elementary  differentials  of 

K  ™~  T-r-Trrn,,.  1  •  * 

order  h  with  coefficient  a  (3.$  : 


V-  »! 


whore  (3*17  and  0  haB  the  following  form 

(1.28)  0=  ^ci1^-1^+...  +  Z_ciH^iM^ 

with 

(1.29)  icjk)?|k)  k=  1 , . . , , h  ; 


i«  1 


i=  1 


(k) 

are  polynomials  in  over  <ft  . 


Proof :  First  we  show  that 

(1.3°)  [f)q  k- 1 , , . . , k  ;  i.1,...,n  ; 

M<  w  the  right,  hand  side  of  (1.2?)  is  proved  by  inserting  (l.3o)  into 
the  left  hand  eidr  of  (1.27)  using  equation  (1.28). 

The  proof  of  (l.3o)  is  by  induction  on  h: 
v-1 :  hero  k= 1  and 


(1)  (lj2on)  i.1  1  _  r  , 

Ri,,  -  lD  yJo  =  <-f  -  * 


induction  from  h-1  to  m: 
let 


(l.3l)  M_]_“  s  _ P  '  *  ^  ^  i,  c  (  l-1»...»|i  5  p»1|...|H-1  )  ) 

’  ^  ord  F—  (j. 


then 


R 


(k)  (j.Pob) 

i ,  H-k 


K-k  . 

»-"■  1  — 


— t! 


x—  1 


(«-lQ  I  ,  r-  flL 

v-v  ai, ayT 

I  x  .  T  u 


Dky 


(1.21) 
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M 


(1,32)  /iVV  e,)  b(«T)\H(«,)  g(«»)  raLDky  1 

Up‘"  «i  J,'VV  Jo 

Inserting  hero  the  induction  hypothesis  iso  obtain, after  some  computatio 
a  linear  combination  of  expressions  of  the  form  Dky|  » 

L 1  T  ayT  Jo 


dy 

where  F^^  is  of  rrdor  x^  . 

In  Lemma  2  (p.®°  )  wo  shall  prove,  that  the  above  expression  is  a  line  * 
combination  of  elementary  differentials  of  order  (h1+...+x  +k)-  x  wher 
the  coefficients  are  natural  numbers.  Done. 

Corollary;  it  holds  that 

<’  *34)  (*].  • 

Proof  jinsorting  (1.28)  into  (1.27)  we  obtain 

v5"  T^0(k)  *».  n(k)  «  y  bT*  y‘o^(k)fFl  - 

fenSr'1  T^7TRi,H.k-  '  1  LJ° 


r*n(k) 


orti  F-  x 


ord  F»x 
,(k)| 


M( 


Hence  follows  (1.34). 


This  theorem  suggests  the  following  definition  of  so-called 
"  elementary  weights”  (Butoher  /  1  /  p . 1 94 ) : 

Definition ;  To  eaoh  elementary  differential  F-  {F^...FgJ  corresponds 

an  olementary  weight  0-  ,  where  0  is  the 

corresponding  coefficient  of  F  in  (1.27). 


Remarks : 

1)  In  addition  equation8 (1 .1o)  and  (1.27)  adjoin  to  eaoh  elementary 
differential  F  a  number  o  and  a  number  (3.  The  correspondence  of 
a,  p,  0  to  F  shall  be  expressed  by  similiar  indices. 

2)  In  contrast  to  F  the  coefficients  a,p,  and  do  not  depend  on  the 
function  f  of  the  differential  equation  (1.1). 
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The  next  task  is  now  tho  computation  of  and  p: 

for  this  wo  insert  (1.34)  into  (1.32)  and  use  the  formula 

/M(x  i“tfi)!  i 


V 

o(k) 

Ri , H-k 


V 


H-k 


v"  v* 

x-  i  h.  +. . . +x  =H-k  ord  F.  =  k.  ord  F  =h 

hJ»  1  T  11  T  T 


t!  x'“*x  ! 

1  T 


(h-t)  1  P1 


/  1  _S  •  (^j)  (  ® .  )  \  /  ^  1  ^  (<j)(d)v_  jT  .. 

VS  yJ< 

I  i  T  IT 


H-k 


(1.35) 


t*»  1  H.+.  .  a  +H  aH-k  0 

hJ  >  1  T 


•  z. 


^Tc  ord  F^=h^  ord  F 


—  ^ hi- 
XT  ^  ' 


with 


[v-%fe*3 


dy 


/  \  M1  1  51  (o.)  (c. ) 

^  ^  ^  (i-1 . t)- 

Tho  insertion  of  (1.35)  in  (1.27)  gives  now  with  x*r: 

* 

y_.  r o!kJ  >  t  y  ..<  _ i_  tt(  i i  t*,?;  4|. 

k=1  f^T  t^T  ri+..  ,+r  »r-k  ord  F^-t,  ord  F  »r  T‘\l»1  K1‘  1  J 

1  x  i  T  11  T  T 

(1.37)  rl* 

.  If  «  dT  _k  o 

•Ifr-F  —  Dyj  . 

dy  o 


Comparing  the  right  and  left  side  of  thin  identity,  we  now  obtain 
recursion  formulas  for  $  and  p:  for,  on  the  loft  side  P  and  0  correspond 
to  F,  whioh  is  of  order  r{  on  the  right  side  p^  and  correspond  to 

Fp  whose  order  is  smaller  or  equal  to  r-1. 
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Conditions  for  tho  Parameters 


-‘or  m-1  in  (1.15,16),  hence  for  the  classical  Runge-Kutta-method  ,  the 
comparison  in  (l.37)  is  easily  done: 
becauso  of  (1.17)  *0  have 


-i,l-  piKl  •»«  “P5 


>o  for  k=2 ,3,...  • 


Hence  (1.37)  becomes: 

WMo  - 


y 
ord  f-r 


_ j _ -x _ a(tf  £«)• 

ET  t*T  r1+...+r  «r-1  ord  F1-r1  -rd  F-r  \1'-1  rl*  j-1  J  ’ 7 

__  I  v  4  H  •  I  H  H 

ri^  1 

.[{ *V”V]o 

If  tie  choose  in  the  left  hand  side  a  fixed  elementary  differential, 

IH  O’ 

say  F- {  F1  . ..Ffl  }  with  j^-s,  this  appears  in  the  right  hand 
side  with  the  following  coefficient : 

-1  US' P"  upyw-K  ■ 


The  factor  "-t * 1  is  due  to  the  permutation  whioh  are  to  bear  in 

Hi  U..MJ 


•*1 

mind  in  tho  sum  ^ _ 

ri+...+r8-r-1 

Comparison  now  gives 

i 


1) 


0.3*)  -.^gr 
o.3*)  tpPt&WT  ■ 

These  are  exaotly  tho  formulas  of  3utchor  /  1 /,  p.  194-196. 

1)  The  coefficients  P  defined  here  differ  from  those  of  Butcher.  His 
coefficients  are  obtained  by  puttin  P’-p/r  • 


v.i. 
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For  m>1  this  comparison  is  more  oomplioatod.  First  we  need  the  following 
two  Lemmas : 

Lemma  1  :  Let  be  functions  of  , . . • ,yn  and  sufficiently  often 
differentiable {then  with  xii+.*»+Hsi  it  holds  that 


(1.59) 


dy,  ...dy,  (W**ut) 
“  s 


-K1 
d  u 


xi - •  •  *  x — . — v 

11  It  al  at  dy,  ... 


1 


“t 

d  u. 


Bu*° 


"al*  0 


9y 


‘si 


i 


-  x1t  .  Mst 
dy,  ...dy, 

J1  Js 


Proof  (by  induction  on  s) : 
s-1  :  a 

w  S...O-  Z— 

h  Hii+... 


j"1 

d  u 


1 


*t 

d  u 


t  ^  (  since 


+(L  ,  «=  1  *  11 

.  .  ..  „  It  oy, 
Hn^.O 


9y 


l1t 


du 


1 


9y<  2" 
J1 


9u, 

li  ¥7  u5' 

j2 


dut 

t  ’  •••  +  ur,,ut-i  dTT  • 


Induction  from  s-1  to  s: 


ds  f  v  d  f  d3*1  /  A 

dy,  ..dy.  1  t'  dy  \dy,  .  ,.dy.  1  t'y 

J1  Jg  L'l  V  Jt 

(  by  induotion  hypothesis  with  kU  m  +...+H  ^  ) 

,K1 


Kt 

a  xu. 


H2i>° 


Msi^° 


yJilxo.+...+H«.«1  V.  . ,+K  .-1  -  K21  a  "Si 

J1\..21v  .  2t  ..si  v  _  st  dy.  . . .  dy  j 

t)** 


a  2t  .  *st 

...ay,  ./ 

"2  Js 


(by  commutation  of  rr-  with  summation  sign  and  by  (1  •  39)  with  s*1  ) 


aX1+H11 
3  u 


1 


v  tii  v  _ 

Kn  +  ...+Kit=1  Mslt,:,+Hstel  dy/1.  ..dy*®1 


>t1i>  o 


Bsi*° 


1 


Bt+B1t 
8  ut 

,  B1t  ,  Bst 

SyJ, 


Putting  now  n^«=  Ki+X<;tt=  h1^+...+h? 


si 


’S  V1  8 

we  obtain  (1.39). 


Pone 


f. 
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For  Lemma  2  wo  need  a  further  symbol,  which  is  now  defined: 

H1  K* 

Dofintion:  Lot  F«  f  -  {F,|...Fp}  with  order  r  and 

p«*  k^+...+hii  ,  and  F»  {  with  order  $  ^ 

be.  elemontary  differentials.  Then  we  define: 

(l.4oa)  f#(Hik}F1f...,Fn)-  {  1)*(ulk}F1,...,FJl):.  {F*11..^**1} 


( 1 .4ob )  (Hik’F1 1  •  •  •  #  ) «  {F,j...i^  }#(n^k{F^  t  •  •  •  *Fn)  :■ 

1  *1  “*F)t  F1  ^*ik  * ,rt^Kik 

where 

1 .7^>  non  negative  intoaor j 
2.)  set  of  indices  for  : 

(1 .400) 


i*  1 ,  •  •  • ,  it 
k- 1  , • •  •  |  ^ 
1-0,1,...,* 


and  for  h„  : 

-  ’  1  *■"'  •  lie 


(l.4od) 


(l.4oe)  k|[U,  with  m-  V  3  ,+k  :»1)  . 


i—  1 ,  • • •  ,  ji 
k*  1 ,  • .  * ,  r 


ik  ira 


j-o 


Convention:  If  misunderstandings  are  not  possible,  we  shall  write  for 
^(Hik,Fi,...,Fn)  shortly  F*(xik)  . 

Example: 

F*  {  F.| 1 .  }  and  $-  { {  f  } f }  :  hence  F^-  {  f)  ,  Fg-  f  | 


xft)  k<2> 


2,  1,  r*  4; 


^(K^)-  {f/1  ••*F11k1}  (by  (l.4oa)  )  , 

X(1)  H0)  K(1)  K(1) 

$*/h0)n  ,p*H  pNl  ,  P  12  \2 

F1<Kik  )m  {F1  ' * *Fn  <  F1  •••Pn  >> 


1)  The  quantities  $,  S',  ?,  ¥±  correspond  to  F  |  rfe,  c^,  $k,  ?k,  ?k  ± 

A  ^ 


correspond  to  $k  • 


I 

I 


.(•) 


«<?>  ,<’> 


V.l. 

x(1>  „(1)  K(1) 


,(2) 
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,(2) 


$*(Hik)-(  Pi11  ...P^1  (P^1  *..Fnn1  tF^2  ...Fwn2>>  (P^1  ...Fnn1  >> 
(with  (1.4oe)  ) 

x««  h  .  X4„  k  „  x. _  x„_  x..  x.. 

.  {  P1 1 1 . .  -P^^  P112...Pn^ 2i  F^3..^ 3>>  (F^4..^  n^))  # 


n*.  .  /S  A  /\  , 

Lemma  2:  Lot  P-  {F^  •••Fn  )  »  t  >  and  P-  lF1...Ft>  , 

then  it holds  that 


Hi  k  .p 

F,’...P  “  —  % 
1  ”  dy* 


(1.41) 


**1 1+*  *  i +*  •  ^ 


(JL 


X1k*° 


with  p-  x^. ..+xn  . 

The  elementary  differentials  9*(xik)  have  the  order  r+9-1  . 


Proof  (by  induct  ion  on  r  ) : 

*«=1 :  Here  P»  f  and  x^-x^  (i»1 , . . . ,n)  • 

From  (1.4oa)  and  (1.9)  it  folli/«s  that  (1.4l)is  satisfied. 


Jcxt  tho  following  summation  rule  is  valid  s 

'<W-  V*-1™ 

(1.42a)  \  flj  ,i+»..  +fli  ( ,1*1 1 . . .  y k) 

J  J  J+1 


then 

T 

Z _ 


(1.42b) 


„(...)  - 


tf1  +  ...+Oi  +X^+. .  •+\k«<j  <J^  +1  +  ,,,+®i  _,‘1 


...n 


(...)  . 


*«  <,ik+1+-+<,t"Ak 


For  simplification  we  define  the  following  sets  of  tupols: 
^ .43)  (J»1  »•••»?!) and 


Induction  from  9-1  to  9; 

'.,'ith  (1*43)  the  induction  hypothesis  reads  as  follows: 


v 
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(1)  »  ,K(t) 


V'-/  d 

p*i  ...fJJk 


x^! 


(1.44) 


Then 


dy 


with  x^+...+H^-  x^  ,  (l»1,...,t)  . 


F1 1  •  •  4^  Pi . . .F  . . .$/ ^-rf )  U ^ 

1  71  dyP  1  *  dypl  1  11  ' 


2 - r*  •  •  xL _ 

<yiO+***+tf1t“1  0pO+‘**+<,p1i: 


/"../pi  _* 


dy 

?°I 

P 


H(0  . 

-  .(•)*, 


/it  /pt  _a^'j  . . 

1  p  KW  ‘I 

dy 


/  °1o  °po  d*  +*  f\  (2) 

•|pi  -,fp 


dy 


x^'+t 


f) 


4 

x 


(D  x(i)  p) 

n  d 


■« 

dy*  ' 


,*<*>  x^x^x  /x<°>  «<•>.<•>* 


t  «  *  P, 


X  d 


dy 


Ttrt 


•  •  «P. 


n  d' 


dy 


) 


(1.44) 


5* 


,  •••T“  y  y  ,  , 

/  y  / . .  •••  /  —  •  •  •  /  m  / 

.,(.)+  +K(t).K  >51  +>)„  '  K(1)  K(1)  K(»)  kT 

'*1  +*..+Hi  “Hi  \  +,*,+Hw  Mr  *1  Kx  K1  Kn 


-  /-X-  v 


jtnJJi»F 


!  •••  H 


j  . A 

X»,  «<*>  \  /x  rtx  aa  x^  .».»*  ' 

H^!  x^i  ,)*  F1(H^-^(H1V,)F1  ••*Pn  git'  f 


W7  • 


(l»y  r«n.x\rnnging  the  sum  signs  and  using  (1.9)  ) 


H ^  +#*#+H4j  —  r*^  a,| 


T  ...X ... T _  T 

i>-»-  $  x®+_+x$»,  *1 

.{F*1  ...»**  ^(h^))...^(k<*))  >. 


(  using  (1.42)  and  (1.4ob)  ) 


Kf*+K»+... «»>.«,  »t^. •  • +*$ -»« ® k$J 


) = 


(  using  (1.4do)  and  x^0^-  x^  )  n 


V* 


•H 


x11+,.,+x1r»x1  Kni'f»  ••+HnrmHn 


(JJi  V",Kjf!)  **(Klk] 


V.l. 
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(1)  :  This  step  is  verified  by  transscribing  into  the  sum-notation 

(of.  (1 .7*9)  )  then  by  using  (1.39)  and  finally  by  transcribing 
back  with  (1 .7*9)  . 

d  1  i  d  , 

(2)  :  This  step  is  verified  by  collecting  together  in  P1  ...P  p 

all  elementary  differentials  whioh  occur  multiple  iso  we  obtain 

M1*  41) 

F,  ...F_”  : 

I  11 


lot  F, 


'k  Fk  "  Fi»  where  Fk  C{F1,...,Fp)  | 

Hi  i 


Then 


*4. 

(1 »45q)  ^  (i-of ...  ft) 


and 


1-1  l 


(1.45b)  «1  - 


(J) 


Finally  we  comprise  the  sums  with  of.  ,  and  obtain 

n1 


v  V" 

^ _ •  •  •  2_ 

d,  +...+0.  x— i 

**-|0  K«|X 

Jjc-|i  ^  o 


T~ 

"  /- 


C,  +• .  .+(J,  . -1 

K  0  K.  i! 

M1  H1 
0.  '  >yO 

1CM11 

Hi  ! 


(4 


(using  (l.45a,b)  ) 


Analoguously  the  exponent  of  Fgf...Fn  and  the  corresponding  sumo 
are  comprised. 


Next  we  prove  by  induction  that 
.  n 


F^^ik)  has  the  order  +?,  (where  Xj»  ^xjk  ); 

induction  hypothesis:  ^(x^)  has  tho  order  ^.x^r^  (  i»1»...»t) 


then 

A. 


H(°)  H<°)  /  V  ,  V 

F*(Mik)»  { F 1 1 1  ...F^*1  F*(Hik  )*«»Ft(Hik  )J  ha8  the  order  (by  0*6b) 


&  P'S  *<-+«■  pp'S, «  (,i45b)  pf}  *  • 

M1t  n 

finally  F-  (F^  ...Fn  }  has  the  order  ^T^r^  +1*  and  hence 


n  t 


7 


>(Hik)  has  the  order  r+r-1. 

With  this  the  proof  of  Lemma  2  is  completed.  Done. 

<3 w  k 

Corollary :  F^...F„  ~])  y  is  a  linear  combination  of  elementary  differ 

'  '  ” . ‘ ' ;  '  1  w  dy K - - 

tials  of  order  r,j+...rx+k.  This  has  been  used  in  oonnootion 

with  (1.33). 

We  are  now  ready  to  prove : 

Theorem :  The  elementary  height  0-  [0.J...0J  corresponding  to 

F-  (F.J...F  }  is  determined  recursively  by  the  following 

formulas : 

0  »  corresponding  to  f ,  is  given  by 

(1 .  46a)  0  -  *>  cP)  I 
i-1 

0  i  oorrsponding  to  F,  is  given  by; 

(k-Di  (Ki\  JK 

,+V-k-T  V’”V  1,1  ’  ‘ 


cr.- 


(1.46b)  0-  [0r..0J  -  21  z>r  T 

v  k-1  i-1  1  L+777 


.17 . t 

^  ^  0 


n  36)  \y(°)e  V  1—  V"  b(0)y(°)  1*1  t 

,1.5b;  *1,1“  2-  0!  ^—Ti j  ¥l,j 


where  }  1 )  oorrspond  to  0^  (of.  (1.28)  )  and 

) 

£t  05  fcTiJ  ‘1,J 

U1 

The  coefficient  p  corresponding  to  F«  {  F^  . . . F^  )is  given  by: 
S—  -j  I P  j  \  H  i 

(1.47)  P=*  r!  1 1  — ,  I—*1,)  d  where  p.  corresponds  to  F.. 

J- 1  0  0 


Proox ; 

1.)  (1.46a)  follows  direotly  from  (1.27)»  since 


.(1)  (l.2on) 

i,o 


f-JD 


k=1  and  hence  k=1  and  R}’1  [fJo:  so  0“  ^Z°i  *  • 


2.)  Let 


**ii  ^10 

F=- {  ,  F^=  (JP^  ...F^  $  t  F^=  (...)  and  so  on,  until  f 


is  reached. 


V.4. 
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Out  of  all  Fik,  Fikl,  ...,f  we  now  oolleot  these  elementary  diffe¬ 
rentials,  which  differ  from  and  denote  them  by  * » • • v . 

Then  each  of  the  elementary  differentials  of  which  F  is  oomposed 

h  k 

can  be  written  in  the  form  { F^ 1 . .  »$nn)  where  o. 


Example : 

F-  { {  f }  <f  {  f 1 h } }  : 

Ff  {f  ),  F2-  {  f  {f2  }>  ,F3={f2>  ,  ■  F4-  f  ; 

then  for  example: 


V  «Fir2F3P4>  • 

;>.)  Let  M(FjF)  be  the  following  set  of  tupels: 

(1.48)  K(?|F)  :*  {  , . . . ,  h^a)  |  F  » •  •  • • 

Choosing  one  tupel  out  of  M(£jF)  and  putting 


(1.49a)  Hi1+...+Hi£,  .  (i*1,...,it)  and  p-  h1+...+kji  , 

we  obtain  from  Lemma  2 


(1.49b) 


F  +  see 


4.)  From  (l,1o)  we  have 


(1.5°)  F,...p .^rikr=  Z___.«(v*-Px  ^f) 

1  *  dy*  or3  v  1  x  iyh  ' 


.)  xi.'.  .uix.nr i n<.u  <jf  (levrnn  in  ( 1 . 57 )  is  possible: 

In  the  left  side  of  (1.37)  we  choose  an  arbitrary  elementary  dif- 

ferential  F=  {  F^  .  ,.F  ^Jof  order  r. 

In  the  right  hand  side  of  (1.37)  we  now  are  interested  only  in 
those  terms,  which  contains  F: 

for  this,  we  insert  (l.5o)  into  (1.37)  and  let  the  summation  run 

only  on  the  tupels  of  M(f*}F)  j  so  we  obtain  by  bearing  in  mind  the 

number  of  permutations  rf  the  sum  and  by  (l.49a,b) 

r^.  •  ,+r^ar-k 


i 
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*jL! 


1  Hd1 !  •  •  •  HJk ! 


JBl 


*  H. !  •  •  •  K  !  *L  Jo 

I  K 


_ JL _ 1  lyi  0.49a) 

K  J1  *  *  *  •’‘jk*  '  ^  Jo 


,(o)\K)A 

J.i, 


Wo 


Putting 


0.51)  Aik>-  ^-SZ_  7|-(|1h-i_  /!i  t(o)\KJi'\ 

Xi  ord  $«k  M(fjF)  lJl  \j-1  Hjl!‘rj!  ^7  ) 

we  obtain 

# Wo-  r'g-  |T^k)A^).  [Fj0  I 
using  (1.28,29)  it  follows  that 

(1.52)  AW.  1 


/if) 

This,  however,  does  not  jet  lead  to  formulas  for  ' ,  since  usually 
the  set  M($;F)  and  therefor  are  not  known.  We  thus  try  to  find 

a  recursive  determination: 

6.)  We  assume  that  the  elementary  weights  01#...,jZl  whioh  correspond 

to  F1f..,,Ffl  of  F»  {F^  ...F^  )  are  known*  then  also  the  followinr 
sets  are  known: 

(i.55)  x(1)(f’1iFji)-  {  ("^...k^)  |  f^x^f, . f„).  ) 

(1-1 . 5  ,  i  i  ) 

and  hence  also 

,  A  (k)  „  f  00 

(1.54)  Ajvi  r  \  fJlfi  * 

Using  (1.55)  it  is  now  possible  to  construct  the  set  M(F}F): 
because  of  (1.4ob)  we  have 

?*(«ik).  (  f""  ...f^°  ?;(*0))...?;(h(‘)))  . 


83 


v.i. 

Since  we  want  P^...P(J<,)«  F  ,  it  follows  that: 

^1  (Kik^)*  Fj  *  ***  *  Ft  ^Mik^=  Fj  where  1  • 

1  ii 

Since  the  tupela  (...,x^,,..)  are  known  from  ( 1 . 53 ) *  the  still 
unknown  numbers  x(?^are  now  obtained  from  the  comparison 

■ft)  h<°V  ,,  . 

(1.55)  {F^  ...F/  Fj  ...Fj  >  =  {F/.^F^J 

as  follows : 

Let  q^o  (i=1 ,2,. . .  ,c)  be  the  frequency  of  i  oocuring  in 
»  then  (1.55)  gives: 

(1.56)  *[°Jm  « 

Fenoe  it  follows  that  only  those  (j.| ,  •  •  • ,  oan  occur  in 
(1.55)  for  which  x|°^>-o  (i~1 » . .  • ,  Jt) 


Thus  the  set  M(?jF)  can  now  be  written 

0  „'o)  H(1) 


1.57) 


m(p8p)-  <  U.  {  (nfr,...,^ 

^  j  1 »  ... »  j^=  1  '■ 


Mi1  ^a,ii“9i  ^  0  (i-1  *  —  *7l)» 


K(t)  )| 

Kik  •" '  I 

(  — w±k^*~  )«MW(^X 
(1-1,..., t) 


/,)  From  (1.13)  and  (1.12)  it  follows  that 

,9 


_  „  d_ _  _k 

Jr “ps  TT  35  y 
dy 

k-1 


T_ 

•  •  • 


^  y _  i.  I  •  \  |  %  j  -  U1 

t=1  A,j+. . .  Xj.=k-i  A1'***At'  ord  ord 

^  >  1 


Jot,... at  . 


(1.5Q) 


*  Fi*”Fs  ~T  (V**Ft) 

dy 


V/e  insert  (1.58)  into  (1.37)  and  let  the  summation  run  only  over 
the  tupela  of  M(FjF)  jdoing  this  we  use  sgain  (1.49)  and  bear  in 
mind,  as  in  5.)  the  number  of  permutations.  Similar  computations 
as  in  5.)  now  lead  to 


64  V.l. 


#M0- 


(by  rearranging) 


here  we  have  used  the  following  symbol  (of,  Grttbner  -  Bofreiter:  Int<r;« 
graltafel  ) 

(1.59a)  (a;-1;X)»  a(a-1 )  ...(a -X+1 )  (\p1,2,...  ) 

(1.59b)  (a$-1jo)=  1 


The  comparison  of  the  ooeffioientc  in  the  above  formulas  gives 


(1.47) 


0-  r!  I 


"  2_  LJU  o 

1  VV' 
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0-  ±  z _ , . 

k-1  i-1  1  t«1  X.+...+X*- k-1  V*V 

ijn  ^ 

^ }  H:l) 


Those  are  the  wanted  recursion  formulas  for  P  and  $» 

9.)  Formula  (1*59)  can  be  simplified  with  the  help  of  the  following 
transscript ion: 


(1 .6oa) 


0-  C0  V?...#] 

✓  /  '  i'f  °S. 

'  )  I  '  V' 

Si  '  N 


'  i 


with 


?«  -  C01-4.0,] 


<2.  4=1 

(1 ,6ob)  s«  2_..4i  »  z_|ik  (i«2,...,o)  ,  it,,:-! 


i«1  *  k-1 

How  formula  (1.59)  becomes 

0.46)  [2,.. .2,]  - ±_  ±4k)H 


k«=1  i«1  H^+.t.+H  -k-1  1  s 

Ml*°  8 


(k-1 )  1  ^*s* 

x1!...x  !  *1 ,i  •***a,i  ‘ 


Proof : 


a)  to  prove  this  v.  ;■  need  the  follwing  two  summation  rules: 

First  Rule: 

Assume  that  out  of  (A, ,  •  •  • ,  X,.)  t  numbers,  say  A.,  ,...,XT  are 
distinct  '.nd  that  A ^  (i-1,...,r)  occur  6^-times  in  (A,,...,Af)< 
if  further  in  (i^,...,i^)  all  i,.  are  distinct,  it  holds  that 


(1.62)  2 


(a^»...»a^)  1 


a,)  <v 

Ai  •  •  •  A, 


1 


6  '...6  ! 


z 


(X,)  (At) 


!  n - T~\  1 

t  (i1f...,lt) 


Aj  ... A.  • 

1  Xt 


Hero  the  symbol  2 _  denotes  summation  over  all  permutations 

(...) 

of  (...)  . 

Second  Rule: 


(...) 


J^»...»j^*1  1  A  j.,  ^  ...  (  i  d  (j^»..»»d^) 


(...) 


(1.62) 
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b)  In  (1.59)  we  now  fix  k  and  t  and  ohoose  a  tupel  (\1 , . . . ,  X^) ,  which 
satisfied  the  same  condition  of  a)  •  Then  it  follows  from  (1.59)  bea¬ 


ring  in  mind  the  permutations  of  tho  sum  21 


X.j  +» » •  +X^.«k-1 


l  (oh^r^:. 


(1.63)  l:-!0  •  *,1.!.#,! 

(x,)  (*t> 

o)  To  show  the  equivalence  of  (I.46)  end  (1.59)  it  suffices  to  fix 
k  and  t  and  to  sum  up  in  (1.46)  over  the  following  tupels 

(M1  f  •  •  •  »  Hg)“  (°  9  IM  9  9  »M  9  0  9  HI  9  9  «H  9  0  9  9  9  Ml  9  O  9  HI  ) 

(1.64)  1  k  t 

with  «  X^  \  1  and  1  <i^<  ...  <i^  <J 

The  resulting  expression  has  to  equalize  (1.65)  : 

With  (1.64)  it  follows  from  (1.46)  : 


0  l~f  Z 


i-1  1  1<i1<i2<»ALt4:s  (X1f...,Xt)  1  *-  t 

.(O  » 


(k-Q»  ?>(0)  }(X1*  $(o) 

X,!...A+!  M,!  — U.i' 


A  X  'kt  /  A' 

."fi.,i'.,?s,i 


(1.61) 


(1.65) 


^„(k)<*r  7"  -JtelU  -  1  o)  ^(Xl) 

/-rci  irr~z — u - 71  V- V  *  a/. ..6,1 

i«1  lii^.  ,.<i^  s  (i1».M*it;  It  1  t  1 


(  l/l j  J“1».».»t) 


a(°)  a(\)  j°) 

Ifl  W  & 

,Yl,i,MWi.,i*"¥s,i 


Next  we  turn  over  to  the  notation  with  by  putting 


Q 


•)  ,  k« 1 , . . . ,  t 


—  ->  where  (1-1, ...,d)  (  of.  (1.6oa,b)) 


In  (i^,...,^)  all  numbers  are  distinct,  but  not  in  (j1»«..»Jt)» 
since  1  i  J1f ...  i  <3  .  Again  denoto  by  ^  (l-1,...,fl)  the  fre¬ 
quency,  with  which  1  appears  in  ( ^ » • . . , j^) • 

Prom  (l.6oa)  we  see: 


V.l. 
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if  q^>o,  then  exists  a  k  so  that 
^  ^  i)^*  *  *  ^i)t+Q  ^^1+1  (l* 1 » •  •  • »  '  and 


Thus 


,  hy  the  substitution  i^  — >  in  (1.65) 


Z. 


1  ^  i.  (  . .  •  <i,  £  8 


(...)  changes  into  IjL _ _ . . .  (j*d)  (. . 

U  ^  i  ...  Vty 


and 

H 


(...)  ohangc-s  into  XI - q1  ! . .  «q0!  (.. . )  • 


(i^  * . .  •  i  i^)  •• 

Finally  with  (£  j)  h !  =  0*65)  becomes 

*P~.(k)  (k-1 ) !  _ 1 _  V~  T~  /  ,  v 

s  C  ,  ,  |  •  I  /, . ..I,.—  I  (  fi«  I  “1  I  Q.  )  * 

.  1X  L1...AJ  0J...6  i  ,  ,  .  .  -  .  .  M  /  .  4  \  1  1 

i=1  1  t  1  T  1  i  ...  i  (J 


JQC) 


L(oJfrQi  L(o)\ V9a  JV 

VI, v  *“V0»V  •fi1.i,,,fJtti  • 


But  this,  using  (1.62),  is  equal  to  (1.63). 
Thus,  the  proof  of  the  theorem  is  completed. 


Done. 


Remarks  : 

1)  From  (1.4^)  it  can  be  seen,  that  the  correspond onoo  between 

ant3  one  to  one. 

2)  For  m*1,  henoe  k»1  and  h1-..,*mq*o  the  formulas  of  Butoher  (cf.(l.; 
are  oentained  as  special  cases  in  (1.46), 


Definition;  A  method  for  the  integration  of  ordinary  differential 

equations  is  said  to  have  order  p  (  or  error-order  p+1). 
if  for  its  approximate  solution  y(x)  holds  that 

jr(r)-J<*)-  0(hp+1) 

for  each  solution  y(x)  of  an  arbitrary  differential 
equation. 


(1.66) 


Theorem:  A  Runge-Kutta-proceac  ht.n  order  p  if  and  only  if  the  oooffi- 
oients  satisfy  tho  conditions 

1« 

(1.67)  0-  —  -  ^  for  all  olomt rotary  differential!  of  order  r$p  • 
The  constants  y  satisfy  the  recursion  fortuula» 

(1.68)  y«  ryj1...^  ^ 

where  v  corresponds  to  0-  f 1  and  y^  to  ^  * 

Proof:  The  assertion  follows  from 


y(x)“  It  21 _ «KL  (of.  0*3)  and  (l.lo)  )  and 

k-o  k!  ord  F-k 

y(x)-  - ftfOL  (cf.  (1.22)  and  (1.27)). 

k-o  KI  ord  P-k  ° 

Thus  tho  Taylor  series  coincide  up  to  order  p  iff  ( 1 .67)  is 
satisfied • 

Formula  (1.68)  follows  from  (1.1 1)  and  (1.47)  : 


Done. 


Examples  of  Conditions 


In  many  oases  formula  (1,46)  for  the  elemtary  weight  oan  be  simplified* 
For  this  we  give  some  example,  which  shall  be  needed  in  tho  next 
sections : 

We  put 

(l.7o)  a±-  2j>iV  5 

J-1  0 

this  is  motivated  from  transsribing  (1.1 5)  to  non  autonomous  systems 
using: 


1 ) 

'  The  coefficients  y  are  also  tabulated  in  Butcher  /  l/» p • 1 91 -1 93 • 
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0.15')  4k)-  0>ky ) <  vaih  ■*0+h  £»{  ]M 1 5 +- •  •  <5? J^M"0 ) 


with 

0-4') 


Kxcmploa: 


(1.71)  t  -  ±o['K  1 

i-1  1 

(1.72)  [/k]  -  £  fLo^>(k,.15x.l)  ak"K+1  .  ^ 


x®  1  i-*1 


(i.73)  it/)- 1  >«)  m?-’- 1^777777 

71.74)  -  ±  . 

m-1  i®1  1  1  \<J-1  Ol 


e+1+2 ) (1+1 


(1.75a)  r01-..MkJ  -  t  £cW]£  M(ki-1,9)a^f  <*-«>. 

K®  1  i«  1  0*  O'  / 


[r+k)y -,..*/ 


with 


(1.75b) 


Proof : 


ai  n 


0-  [^...0  J  »  T_  £'c[>l)'i,[H)  and  y-  ry1...yfl  . 

w*  1  i= 1 


1. )  (1.71)  is  already  proofed  (  of.  (1.46a)  )• 

Ahe  quantities  ^  which  correspond  to  /  are: 

(1,76)  *»  ^  ^  1  5  o  (k«2 1 5 1  *  * .)  j  ip  ^  • 

2. )  Because  of  (1.76)  for  ourimat ion  in  (1.46)  runs  only  over  the 

following  tupcla  und  their  permutations : 


v.1. 


vt-1  k-x-1 


(h-Di.J-"*1 


3.)  Because  of  (1*76)  for  the  summation  in  (1,46)  ia  only  over 

tho  tupela 

( X^ , .  • . ,  X  )■  (<J,  1  ,  i  •  •  1 1  ,  Of  •  •  .  ,  o)  with  s-k+1  and  o  <  tf  4  h*1 

and  their  permutations: 
henoe 

m  -  ±  (s  ■* 

(V)  (ki-1  ,x-o-1) 

■  £  £*!•’£  W  . . W  • 

U)  In  (1.75)  we  put  ft*1]  .  Thus  with  (1.72): 

(l|-1|H-1)aJ’x+1  (x-1,2,...  ) 

*M*&Sb^)(1,-1,e-1)a5' 


L-O+1 


>1 

inserting  this  into 


we  get  (1.74)  • 

5.)  Prom  (1.46)  **  have 


(1.77)  ^K‘<5)-  2 _ _  (k-o-I);  . 

1  ».  .*  u  4  !  ^  »  111  S.i 

A.+. ..+X  »X-o-l  ’  S 

a 


X,  £  o 


PcoauBP  of  (1.76)  in  the  case  of  [ ^ •  .^^3 the  summation  in  (I.46 


is  to  extend  over  the  following  tupcls 


...  a  k-0 

with 

J- 

X. o  (i»1,...,s)  ,  2— ^  • 

1  i-1  1 

Then 


p,  •••*/]  - 1  £4H)t 


La 


H*=  1  i*  1  1  X..+*  •  .+X  -tf*®  H-1  / 1  1  \  X., !  • 

Xjio  ,  <si  o  S 


<X1>  T(Xs)ak-B 

•  *1,1  *1 


tt 


*-  r  /k\  (k-1)!  k-0 


n  n 


X^  ^  o 

(k)^'  (k-1 


(h-0-1)! 


X^  V  o 

<x1>  (x»>  (1.77) 

*  *  1 ,  i  ••*Ys,i 


^  2S£X) 

h-i  r-i  1  o^ 
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V.2.  Implicit  RungQ-Kutta-ProoQssos  wi. 


Nodes 


Introduction 


Ac  with  clasaical  Rungo-Kutta  methods,  also  here  the  following  distinc¬ 
tions  aro  useful: 

explicit  method:  b^'=o  for  j*i,i+1  5  k*1,2,...,m 

semiexplioit  method:  b^y«o  for  j-i+1,...,n  ;  k-1,2,..,,m 
implicit  method :  otherwise  • 

(k) 

In  the  first  case  the  values  g£  '  oan  be  evaluated  recursively, 
otherwise  they  are  determined  by  implioite  equations  which  may 
co  solved  by  iterations. 


Thu  following  thoorem  about  implicite  Runge-Kutta-prooesses  is  due  to 
~utoher  / 2  /  : 


Theorem:  Each  quadrature  formula  (with  single  nodes) 


n 


(2.1)  y(x)  =  yo+h/_otf  (xo+aih) 

i=  1 


oan  be  extended  to  a  implioite  Runge-Kutta-proooss  with  the 
same  order: 


n 


v!>z,olgj 

i=  i 

whero  a^  ,  c^  arc_ 


by  the  equations 

f.  k-1 

pf  iJ  J  “  k 


with  g^=  f( 
the  values  of 


X0+aih»y0+h^Tb 

j-1 


(2.1)  and  b 


id 


id8d)  ’ 

are  determined 


k= 1 , . . . ,n  j  i«= 1 , . . . ,n  . 


Here  we  arc  showing  that  an  analogous  theorem  is  also  valid  for  quadra¬ 
ture  formulas  with  multiple  nodes. 


9  3 
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Quadrature  Formulas  with  multiple  Nodes 


The  following  generalization  of  (2.1) 

(2.2)  y(x)=  yQ+h  o^f(x0+aih)+h2T~ o^f « (x^-t-a^)-*-. . . 

i» 1  i«1 

i»1 

is  called  a  quadrature  formulas  with  multiple  nodes.  Here  not  only 
the  values  of  f(x  +ajh),  but  also  derivatives  of  it  are  evaluated. 

Such  formulas  (  or  special  cases  )  have  boon  investigated  by  a  number 
of  authors,  e.g.  D.D.Stancu, A. H. Stroud  (/4&/,/$o/)  ,S.Filippi  /1 5 /  •••• 
Here  we  restrict  ourselfes  to  the  formulas  with  multiple  Gaussian 
nodes  given  in  Stroud  -  Stanou  /5o /.  Those  reaoh,  similiar  to  classical 
Gaussian  formulas,  the  highest  possible  order. 


The  following  theorem  is  proved  in  /5o/: 


Theorem ;  If  m  is  odd,  the  coefficients  can  be  determined  so  that 
formula  (2.2)  reaches  order  (m+1 )n  ,  where  the  ooefficionts 
.00  . 


(2.3) 


(k=1,...,m)  are  given  by : 

/_  T  cpOfi-lj-ijic-l)*}-1  =  1  1-1,..., (m+1  )n  . 
k-1 


m 


These  coefficients  are  tabulated  to  2c  D  in  /5°/  for  m»3»5  and 
iv ^2  (1)7  (  for  the  interval  )  . 


Implicite  Runge-Kutta-Procoss  with  multiple  Nodes 


The  idea  of  the  following  proof  is  due  to  Butcher  /  2/,  '^he  verification 
of  the  single  stops,  however,  here  is  very  muchmoro  complicated. 

We  first  prove  the  following  formula: 


9/i 
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(2.4) 


(r^^.+r  |-1  in)-  21 _ f  ■Hl  ■  .■  (r1|-1|H1)...(r  «-1|H 

8  h1+,.,+h-hk1i,,,hsi  1  1  8 

.  1  *  _  S 


«i  >,  o 


where 


(2.4a)  (rj-ljn)-  r(r-1)...(r-ii+l)  (n»1,2,...)  |  (rj-1>o)-  1 

I’roof  (by  induction  on  \i)  : 

the  case  n»o  is  oorroot. 

Induction  from  ji  to  n+1 : 

(i\,+...+rai-1  ju+1)  -  (r1+...+ra5-1}n).(r1+...+rg-n)  • 

(r1l-n*1)...(r8i-1i*8).t(r1-*1)+...*(v«6)] 

\  >  O 


2 _ , 


7Z  <ri  iK,  >  •  •  •  (ra  -‘V  (ri-Ki>  ■ 


i-1  H„  +  ...+H  «U  1  S 


*1 +  *  *  *+H  “4 
1  0B 


rB- 


>  r _ 

ia  1  Hj+t • • ^ 


^+. • • +(x^+1 )+. . .+Ka*H+1 


X^  0 


£n_ 

in  h,+... 


+H  -H+1  "1 


...(ra.-1«x„) 


x^  o 


Jil 


/JL. 


Hi) (r1}-1;H1)...(rs;-1jHs)  - 


H.+...+H  -d+ 1  hi!«**hs!  h-i  y 

H1>0  8 

'y  — (h+U1-  (t  ..1  jH4)...(r„}-1  jh„)  .  Done. 

- - r - -?x,!...x  !  1  1  8  3 

H1+.  .  .+H_=(i+1  1  S 

*  v  3 

«!>  o 

Consequence ; 

.  '  3  —  1  .■  \ 

(2.5a)  (k+1-1 1-1 1«>-  (k-lj-ljs)  +  1 21  (  JJ  (1-1 1-1  |x)(k-1  j-1  |8-h-1  ) 

X=0  *>  l 
8-1  .  v 

(2.5b)  =  (1-1  s-1  js)  +  k3I  |  (k— 1 1  —  1  |h)(1— 1 }— 1 } s— H-1 ) 

X=  0  '  < 


V.2. 


(2.5a)  follows  from  (2.4)  for  r^=k-1  ,  r,j»l  >  4°  s  }  (2,5b)  is  (2. 5a) 
with  k  and  1  interchanged. 


For  the  rest  of  this  section  we  assune  the  following  conditions: 
* 

a.  /  a.  for  i  /  k  and  a.  /  o  (i=1,...,n) 

(2.6)  (k) 

o  ic  1 )  •  i  •  ^ n  f  kR  1 ,  •  • .  ,ra  • 

Using  (1,72)  f  0»75)  wo  now  define  the  following  symbols: 
Defintion: 

A(  §)<•—>  y0  =  1  for  all  elementary  weight  of  order  r 


H\)  «->  ^ _ e±V’; C1^-1 5-1  JH-1  )aJ“K  -  £  k«  1 , . . . ,  f 

1  i=  1 

_  _  k 

r (1 )  /  >  (k-1  ;x-1 )  r—  .  (h)  k-?i  _i  i=1,,,.»n 

M!  k  kBl . $ 

*  k=1  i=1  1  1  1,1 

-  t_  •i*) 

v  k  51=  0+1  J  V0/  J 


» •  •  • »  f 


/  (°)/<  k \  m 

«'  I  il’.-fjl  -  z_ 

v  k  51=0+1 


0=  1 1 , , ,  y  m 
j=1 ,,,  ,,n 
k=  1  ♦ , , , » | 


f-2  wi  i  \4t  0! 

5t=  1  i=  1  y p  =  1 


3  =1 ,  .  . . #n 


k=  1  y  •  ,  •  y  | 


Theorem: 


(2.7)  B(|+1J)  ,  C(*f)  =»  E(|,rj)  * 
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Proof: 


5  ‘■-,'*l5  C;J  <«-'  I-.  ,■-«-»(■-,  ,-„«)! 1 

k-1  i»1 


Theorem: 


(2.8)  B(  |  +  m.n)  ,  E(  |  »  nun)  —=>  D(^) 


Prcof ; 

1  1/11  \b(£+7) 

lTi+kT  "  k  \T  *  Uk)  u 


r  H  21o^H)a^K  ((1-1  ,-1  jK-1  )-(k+l-1 5-1 JK-1  )aj)  (2i5b) 

k  k-1  pi  J  j  \  V 

£  . . 

(<j  — >  x-cr-2)  .  k 

m  n  /  \  ,  „/  1 -a . 

XI  Zo  KM  - 

Si  £t  •*  J  v  * 

»«•>  To)  *  -  $ 

3  it  XT°jM^aj’M((l"1;“l5M'l)  “T*  +  0+ij  (k-1;-1jH-ff-2)(l-1}-1j.V. 

-  II  ^cWaf‘,t  (k-1  ;-1 1 H— 0— 2 ) (l— 1  }-1  }<j)  - 

i^T  JTi  1  J  *+1/ 


1 -a*  .  h-2  a 

i)  "jf  “  aj  21 M  (k-i i-i ia)(l-i i-i ik-o-g; 


^Ilh1)  (k’1  J-1  JH-<7-2)(l-1  5-1  r 


( 

\ 


in  the  last  expression  wc  now  substitute  fl  —4  h-1  and  h  — >  (J 

I 

J 


the  summation  bounds  o  (  U  ,<  h-2  }  j-2  thus  transform  into 
Vo  £  m-1  £  <J-2  ^  m-2  ,  hence  K<hJ  n-1  ,  H+1  £  <j  £ n  • 


4 


'*’  "’-“-’X1-1 1*1  I*'1)  * 


m  m+1 


/  Z. _ 

X=  1  0»  X+l 


(  by  rearranging  ) 


m  n 
-c*  v 

•1-i  jTi 


(1-1 }-1 }h-1 )a 


r(»SK)  ^  D«r*} 


(2.9)  m  n  /  \  k-2  .  v. 


h*  1  j-1 


0*0 


(k-1 1-1 >h-<J-2)(1-1 j-1 jo)  , 


on  the  other  side  ee  have  from  E(^,ir|)  after  slight  modifications: 

nkr E(J,,)  t  ga-’.-i.-D-r-  *(£ 

M  (k-1 1-1  |H-0-2)(l-1 1-1 1®)  . 


(2.1c) 


m  n 


*  2_  >  ‘-s: 

x*  1  i*T  1  1  c-o 


Subtracting  (2.1o)  and  (2.9)  we  obtain 


*(° 


(») 


1  -a 
it 


m 


0«H+1 


1<a)C)(k-..-lia-1)a$+«-)).1 


Nov;  v;e  consider  the  expressions  in  the  saved  brackets  {•..)  &s  indepc-nii  . 
variables  and  let  1  run  from  1  to  ra.n  ,  then  this  is  &  homogeneous  linear 
system  of  equations  with  non-  zero  determinant  (  cf.  footnote  on  p.icj 
and  (2,6)  ),  Honoe  the  system  has  only  the  zero-solution  {...}■  o 
(k»1,...,m  5  j«1,...,n  ;  k-1,...,^)  ,  this  means  that  D(  | )  is  satisfied 
Done. 


Defintion :  Given  Z 

H*  1 


n 


.(h)v 


i=  1 


(h) 

M 


then  we  write 


90 


V.  2 


(1.40)  shows  that  JZ^*  ^  (i-1,...,s)  yields  [0<|*..0a3=.[$1...08]  • 

Lemma  1  :  If  C(tj)  holds  and  0  -  [^...0  ]  is  suoh  that  ^  have  orders 
ri^  ^  (i"1 » ♦ • • » s) ,  lihon 

(2.12)  C^’1]  . 

?:*c.pf :  (ty  induction  on  r  ): 

for  r--2  (2.12)  is  satisfied  ,  since  is  the  only  elementary 

?eieht  of  lhat  order. 

First  C(r|)  gives  since  r^Tji  (i«1,...,s)» 

(2.13)  ?  2  [Of11"1]...!:/8  ’]]  5  — ^  [/r_1]  with  «8+1  j 

1  *  s 


because : 


a  fk)  r  .*  i- 

the  coefficients  '  belonging  to  \_p  J  read  as  follows: 


r,-1 


(2. -14)  -  (rt-1 1-1  ik-1  )a3 


r^k  (r^-ljk)  r.j-k 


a j  (k— o ,  1 , . . . ,m) 


^or  k-1,...,m  this  follows  immediately  from  (1.72)  and  for  k*o  fron 
y(o)  .  Y  1  .  Y  -  T"D^H^(r  ,-1|X-1)ari’’<  cJC‘T?  >  “i1. 

‘•i  h  fcr*  ^  h  «'  (  ,x  ,  k  ^ 


Thus 


d  (1.46)  f-  f-0(H)ri _ _  (h-1  ) !  ft^l*  (2J4 

h^t  j  V*,;+Vx’1  *”  s,d 

0^  *  0 

#-  ^„(»)v  (h-1  )  i  <ri»-,»ai)  <rsi-1in,) 

:£t  1..1  J  B.+...+0  -H-1  <V‘*,8s1'  r1  rs 


'iv - s 

*  o 


r-H 


.a  . 
J 


r1+...+V(x-1)  (2>1J) 


H-1  J-1 

(  udinc;  (2.4)) 


1"  s  "IE - s 

0.  ^  o 

i 


Z  .■■■.'p— H-1  (r,  I-1  *<*,)..'  (rfll-1  l[ 


V.2 


99 


1  xr-  (h)/  i  i  r-x  (1.72)  1  r/r-li 

- T-  2_  /(r-1  j-1  jk-1  )a.  v  4'  - - -  \J>  ]  • 

r**ra  h-1  J  r1*'*V  J 


n  n 
'  ^ 
d 
$ 

Induction  from  r-1  to  r  : 


Induction  hyp  0  thesis: 

r  .-1 

(2,15)  Vi0i-riC^  3  »l*ere  ri^  (i-1»...»a)  1 


thus  we  have 


r,-1 


Done. 


Corollary:  If  C(t|)  is  satisfied,  then  (2.12)  is  valid  for  all  elemen¬ 
tary  weights  of  order  r£»|+1  . 

Since  for  those  the  conditions  of  Lemma  are  satisfied. 


Lemma  2 :  If  C ( TJ )  is  valid  and  <f)  =  [01»**08]  is  suoh  that  for  the 

corresponding  orders  it  holds  that  r^if+1  ,  r^lj^  (i»2,...,s) 


then 


r-r.-l 


(2.16)  y  0  5  ry1  [f^  1  ]  . 

Proof : 

First  wo  have  by  Lemma  1  and  C('(|) 

(2.17)  W1[/?'1]...[/r°‘1]l5-J-r  Vj 

2  s 

where  ?  r.^  (i=2,...,s)  ; 


r,+. ,,+r 

2  S  ' 


since 


A  .  s 
a,  ^  o 


(  using  (2.14)  ) 
~  n  .(h) 


y~A  ^  ^~7rar=“Sriol!**-0s!  ^ 

^  «  V  s 


(>.-1)1  (°1>  (rgl-1|fl2)...(r,  1-1)0,) 


r2*  *  *rs 


Iti  o 


v.p. 


. 1--  f-  f  0(h)^1  — —  y  irv2! 

V-r.  jti  J  2^  «,<<»-«; -1)1  .wri  V-V  ‘ 


(  using  (2.4)  ) 


•  (r2l-t|C2)...(r8}-1|<jg)aj* 


r2+. .  .+ra-(K-<J1-1 ) 


r-V  t  fol  (r2+...«8i-1,K.ai.1)a^- 

r2  *rs  M-i  3  p^ToW  28  1  j 

(  <J1  — >  H-^-1  ) 

-  m  n  /  vx-l  .  A  r„+.  t.+r  -0 

r2**,rs  H-1  pi  J  ^  W  2  8  •* 


r2+...+rg-(H-<J1-1)  (i 


r2+...+rB,-1,ff)a  ‘ 


s"%(x-(J.l)  (1J3) 

*  4  j  — 


7777“  Wi/ 


V— +rs] 


2 - 8 

Thus 


«*>«  *V"r.  DV-/*?  (2^2)  *.*•••*,  [Z2'’]... [/’"]] (2^ 


n,  DV 


r2+...+ro 


i/one  • 


Theorem: 


B($)  ,  C(ij)  ,  D(J)  } 
( 9 . 1 B )  f|^+Tj+1  ,  j  $2  17  +2 


j.e . ,  the.  Runge-Kutta-process  has  the  order  ^ 


Proof: 


1. )  If  0=  (>,... j0g]  with  r^f (i*1,...,s)  and  r£^  we  have  by  Lemma  1 

and  B  (  ^  ) 

rX  <2iJ2>  rf/r-’]  ,  . 

In  particular  this  is  the  case  if  r<f|+1  (  Corollary  to  Lemma  1  ). 

2. )  Thus  it  sufficients  to  oonsider  the  elementary  weights 

<p  =  £0.j***0a]  for  which  at  least  one  0^  ,  say  0^  ,  has  order 


V.2. 


1  ol 


r,>Tl- 

From  the  oonditionj  <2i^+2  it  now  follows  that  the  others  elementary 
weights  have  orders  r^  TJ  (i«2,...,s),  otherwise,  if  for  example 

r2^/l‘thien  we  have 

^r-r1+r2+...+rs+1  ^  (v|- l)+(t| 2ij+}>$  , 
hence  a  oontradiotion. 


Thus  the  remaining  elementary  weights  ^  satisfies  the  conditions  of 
Lemma  2,  and  it  is  suff^ent  to  show  that 

(2.19)  1  J  =  I 

since  then  we  have 

4  (V6)  rriGMr‘r,‘]  (2-19)  1  • 

Proof  of  (2.19): 


This  proof  is  by  double  induction  over  r  and  r^  for  r» +2 , . . • 

r1=7j+1,...,r-1. 

First  we  have  from  the  conditions  (2.18)  and  beoauso 
(2.2o)  r-r.,4  | +1-(t^ +1)  «=  ^  5 

and  hence 

<y»ir.wL,, iI,-1)ap-\w. 

h®  1  i*  1  V  * 


and 


A*2C-V.  t  r-r,-n-1  (h-o-1)^  (K}6) 


(J®  0 
n  n 


M 

r-r.-H  m 


ra  n  /  \  /  r-r.-H  m  .  n  /„\  /  \ 

. .  £ ^ 

. • 


H=  1 

h-2 


C' 


(interchange  of  sequence  of  summation  in  the  first  expression, 
substitution  0  — >  h-o-1  ,  i — >  k  in  the  second  ) 


v  -ffl  V'  v-0M(r_r  ’  bWlf(ff)  * 

Si  i&tV”  h  h  1  (  1  *  ’  ;  1 

y  y  c(*)y  |w“1|(r-r  -1  --1  .H-o-1)n.:r"ri"’<+t't(Kl  i (2.2o)  and  I>( 

£ rtf  *  £rWl  1  ’  '  *  '•* 


1  o2 


V.2. 


JL  JL  /o^(l-a.  ^  ® 


K=0+ 


c,  env  '•>,  ]><:, 

£  •  (-« oo  pi 

^ „(«>.(*)  . 

r-r  ^ ^ T1  k 

T  r1  o»1  k»  1  K  ,ic 

4  m  n  /  N  /  r-r4  (k)  x-1  /  A  r-r.-n+o  ( 

(  a  — >  h-0  ) 

i  <  m  n  /  \ H“1  .  jl  r-r.-a  (k-o)  /.  -c\ 

^(V,  -c*,~  •>/"’)) 

with  • 

Since  r  »e  have  the  induction  hypothesis1 ^  and  since  the 

highest  order  of  $^,,,,,0^  is  smaller  than  r1 ,  we  may  also  assume  the 
induotion  hypothesis  for  r^  j  this  yields 
ri  1  yr“r1l  1  r1 

Wr-V  J  “  r7r"..^t  "  * 


This  gives 
r-r4-1 


tty  1  1  ■  7^7(7,  -  577)  -  777  > hence  (2-19)-  Bo”e- 


Theorem:  it  holds  that 

(2.22)  B(m.n+n)  ,  C(m.n)  *»>  A(m.n+n)  • 


Proof : 

^  =n  ,  >|  =  m.n  ,  ^  =  m.n+n 

(2  7  } 

B(m.n+n)  ,  C(m,n)  -  *  -  ■>  E(n,m.n)  ; 


1  ^The  induction  start  with  r=r|+2,  r^'rj+l  where  (2.19)  is  correct 
since  then  the  elementary  weights  in  (2.21a)  satisfy  Lemma  1. 


lo3 


V.2. 

B(n+m.n)  f  E(n,m.n)  D(n)  j 

since  and  2^+2  are  valid  we  have 

B(n+m.n),  C(m.n)  ,  B(  n)  A(n+m.n)  •  Bone* 

We  are  non  ready  to  formulate  our  main  theorem: 


Theorem:  The  quadrature  formula  with  multiple  Gaussian  nodes  (2.2) 
pan  bo  extended  to  an  implicit  Runge-Kutta-prooess  with 
multiple  nodes, This  has  the  aame  order  than  the  quadrature 


formula. 


M 


The  ooeffioionta  are  uniquely  determined  by  the  nodes 


ai  and  the  weights 


00 


Proof: 


B(m.n+n)  is  satisfied  because  of  (2.3)* 

By  _  _  nk 

c («.«)<=»£  ifclisf ifcli  i 


k« 1 | • . t  fDiU 

i-1,...,n 


(k) 

the  ooeffioients  b^7  are  uniquely  determined  sinoe  the  corresponding 
1 ) 

determinant  '  does  not  vanish. 

Finally  from  (2.22)  we  have  that  the  Runge-Kutta-process 

has  order  (pv.n+n) .  Bone* 


We  still  remark  that  all  other  theorems  of  Butcher  can  be  generalized. 
We  state  them  without  proof  since  we  do  not  need  them: 

Theorem:  B(n+  ^)  ,  E(n ,T|)  =>C(vj)  $ 

Theorem:  B(  ^  +  ^ )  ,  D(  ^  )  ==>E(  ^  ,  lj)  )  • 


^this  is  a  so-called  "confluent"  Vandermonde  determinant  whioh 
is  regular  if  and  only  if  the  nodes  are  all  different 
(  cf.  Gautschi  / 1 7/ ) •  But  this  is  assured  by  (2.6). 


1o4 


The  iterative  Computation  of  the 


We  non  show  that  the  values  g^)  which  are  determined  by  the  implicit 
system  of  functions 


(5.23)  g|k) 


CBky)(*o«1h,yo+h|.t(1j) 


,(1)+ 

Sj  +,“ 


( k*  1 1 . . .  ,m  |  i=  1  y  •  •  •  f n  ) 


can  be  computed  iteratively: 
ire  put 

=  max  |b|^|^|  (k«1,...,m)  and 

it  v  ||  =  max  { ,  •  •  •  »vn|  with  v=  (v1,...,vn)  (  vootor  norm). 

k 

Theorem:  If  the  functions  D  .v  (k=>1 . m)  satisfy  a  Lipsohitz-conditicr 

(2.24)  |(Dky)(z»)  -  (Dky)(z' *)  IN  Lkiiz'-z"  il  (k-1,...,m) 

In  somo  domain  B.  and  if  the  step  aize  h  satisfies  the 
following  conditions 


(2.25a)  |h|LB.|0  where  L=L^+...+L 

Ihf-1 


(2.25b) 


k! 


^  B^  (k=2 1 . . .  fin) 


then  (2.23)  possesos  a  unique  solution. 


(k)  -(k) 

rroof :  Assume  that  there  exist  two  solutions  g^  '  and  g^  . 


then 


ii«f°-*ik)i  (2<‘°  iJhtbi1j)(41)-®31))+-^r  i' 

J=  I  d” 1 

Iv(|h|B,max||gj1,-g<l)||+...^£  BBma*||«jnl)-gW||j  (2-s25b) 


^  .1 

Since  this  is  valid  for  all  i  9  v?e  have 

max  ||gp^-g?k)  ||  4  L.|h|B1  Xjnax  !|  g^-g^H  and 
i  1  11  K  1  1-1  .1  J 


V.2. 


1o5 


m 

X 

k=1  i 


max  ||gjk^-g{k^  ||  i  L|hlB1  T" max  ||g^-g!1^||  with  L-L.+.  ..+L  . 

j  A  1  j  J  J  *  ® 


This,  however,  gives  a  contradiction  with  (2.25a)  , 


Theorem;  Under  the  conditions  of  the  preceding  Theorem  the  following 
iteration  converges  to  the  solution  of  (2.23): 

jlL  *-„(«)» 

j’ 

(k-  1 , .  •  •  ,m) 

with 


(2.26a)  gW.  (Dky)(xota1h,y0+h|_b^)g(^.1+.... 


■  !  13  ®3|B-1; 


(2.26b)  g<k>  .  O  ,  -  [nky]0  . 


Proof; 

With 

(2.27)  K  =  L|h|B1  < 1 
wo  obtain  analoguou3ly 


H«Mt  "eMI-1  II  4  LitlhlBi  Ilmax  ~g J^N-2  H  and 

(2.2S)  Xlmax  ||g{k>  ^max  II  B^n-1  “*i^N-2  0  (N-2.3,.-.) 

kc  •  i  k*  1  i  ^  ^ 


Thus  we  have 


« «$  i»i>,  p»  n(2i27) 

^Kl?max||gW.1.6(l).2||  (2^8) 


Lk 

T 


K  fTi””  ^gJ»N-2“ej,H-3  ^  * 


(2.26b) 


L,  ,T  .  m  , 
U  1-1 


const  K 


IT-1 


||_(k)  (k)|.  c  ||  (k)_  (k)  .i|  ,|  (k)  (k)  j|  . 

l|Si.N  gi  11  *  I' gi,N  gi,U+1  >1  l|gi,N+1  gi,N+2'l  *** 


1  ;>C 


V.?. 


const  1-K 


Done. 


Table  of  ooeffioionta  for  m»3  »  n«2,3»4. 


The  nodes  and  weights  ojk^  are  those  of  /48/  transformed  to  tho 
intervall  C°»0» 


The  coefficient s  are  tabulated  in 
the  by standing  sequence. 

Condition  (2.23b)  gives  for  tho 
step  size  h  the  following  restric¬ 
tions  : 


n-2  : 

h  <11 ,9 

n«3  : 

h<15,3 

n»4  : 

h  <2o,0 

This,  however,  is  no  limitation 

to  tho  practical  use  of  these  formulas. 


Order  8 


tlP539*435325o45/foo 

tni46o55Ml7A95i/*oo 


,300000000000000/ 400 

,24o729*2o844974/-ol 
,  36626496o671727/-o2 


,2olfl5*115ftHoo5/»oo 

,5l^596(?ooo5959/»oo 


223466569o8o54 l/-ol 
,5683467l8998l9o/-ol 


,U67396684oo997/-ol 

,24i294lol5o96l5/*ol 


,500000000000000/ *00 
-,24o72942o84^<i74/-ol 
, 36626496067 1727/-o2 


- , 16^59 6S00059598/-0I 
,298l45^4i63994/.oo 


, 868878778082417 / -o2 
-,7o49254lo77o49o/-ol 


- , 21535125I065784/-02 

, 1o3o193o8oo2o39/ *ol 


V.  2 


1o7 


,927%4o72111lR3/-ol 
,  300000000000000/400 
,9o7?.195927flM8l/4oo 


,2^)6582029^!'’/ 400 

,7791l6664923388/-o2 

,313435o</115744o/-o3 


,1<'3773(,H343513o/40o 
,  37  39352066093  /  400 
,  2.627  93Mo7  1^598  /  400 


-,4'6591o9667o57'V-o2 
,  179  65259754  5^  25/ -ol 
,  1  495R37  2  3.V  2^  7  / -ol 


,  l65925i3959o37l/ -o2 
,  333750^  5oR69p6/-o2 
1 294 55n526ool9/ 3/-o2 


,55iui593l9722o5/-ol 
.32o37o736>' 80558/ 4oo 
j  679<!  29263319'!4l/4oo 
:944838'iort8o2779/40o 


,1/234?/W  8809491/400 
,3oo430;)c<)36l.899/-o2 
fil2't67(3l9565°7o/-o3 


,623726379874  444/  -ol 
t\CT!^mHr)7%\l*no 
,•6005*946718123/400 
,l636377796396o6/4oo 


-,14o97264 5945651 /-02 
, 68727° ’95610269/ -o2 
, 562012618833977/. o2 
,624ol7o67?.33256/-o2 


,365378o96>l454/-o3 

. 733089; 33578853/-03 

,  64693  ’7'>9'8A99^o/-o3 
,691 49749068'' 649/-o3 


Order  12 


,466683594o78222/4oo 
, onoooonoonooooo/ 400 
, 276598562227198/-o2 


-,1486824o47o3o24/-o1 
, 23334 1797039111/ 400 
,481551934548324/400 


,275196431 284 215/-02 
-,387o49R437  !94:W-ol 
, 275196431284215/ -o2 


-,1992374o'!6llol8/-o2 
,92976568663l594/-o2 
, 185876877707 l2o/ -ol 


,266658202960888/400 

-,7791l666492838R/-o2 

,513435o9115744o/-o3 


, 39?476224629o8l /-o2 
-,fi38ll491o32o6l4/-o2 
,l62R3431752575n/*oo 


- ,  lo96496o651^!297/ -o2 
, 238292645997474/  .02 
- ,  2oo432'!  426527 35/-ol 


,1351o52R6925158/-o3 
*  1 276897903924422/ -o3 
,1421359151o4o92/-o2 


Order  16 


,337656551190308/400 
,234766n<>2367929/ -o2 
,  lol669;!  56137282  /  -o2 


-,llo825972275662/-ol 
,l6i;)576>’2l47363/40o 
,  34  3':-74lolHo6443/  400 
,33?.^o667923o5l/400 


,7o22563R4227°46/-o3 
-,1717719929o6364 / -ol 
, 69629355574 o32o/-o2 
,  129126839112262/ -o2 


-,774766653o59'»72/-o3 

,31o534o988o966l/-o2 

,  6781 733  39?754 8/02 
,  574o76l7‘i  33o44  2/ -o2 


,33765655119o5o8/400 

- ,254766802367929/ -o2 
,  lol669456l37282/ -o2 


,  521588326745738  /-o2 
- , Io8l755o6l59353/ -ol 
,1762989o9o42644/400 
,34n739l4n4l8o74/400 


- , Po4 067656235976/ -o3 
,1?6739951oo4463/-o2 
-,222725353379949/-ol 
- , 4  393o791 631 3o57 / -o2 


,3594o56244  32543/ -o3 
-,6Bl566o7!133556/-o3 
,299'32639o14o32/-o2 
,  687‘  934o2129687 / -o2 


,162^34488o9491/400 
-  ,3oo4382o936l899/ *o2 
,H24678l956587o/ -o3 


- ,  134433o83o115o5/  -o2 
,22925o2o9136777/ -o2 
-  ,5394694o47769fl8/  -o2 
,9997o8lo822o474/ -ol 


,2314o6485o94629/-o3 
-,  338637998998192/ -03 
,86402076886*1724/ -o3 
-  ,74l349o6466945o/ -o2 


-,  1689o573289429o/ -0 
,  2787218754523^17 /-o 
-  ,582825o6l3363ot/  -0‘ 
,3o942882o59377l/-o3 


V.3.  Explioit  Prooesa  of  Orders  m+s  (a  $  5) 


In  this  section  we  shall  give  some  explioit  prooess  with  multiple  nodes. 
General  theorems  on  their  existence,  as  with  implicit  methods,  are  not 
known,  of  course.  One  has  rather  to  content  one  self  with  a  laborious 
search  for  special  solutions.  The  question  for  "  optimal  "  methods 
(few  nodes,  small  error)  is  still  more  diffioult. 

Thus  wo  assume 

(3.1  a)  a1  -  o 

(3.1b)  ^ij^-0  }  k*1,...,m  )  . 

"’o  further  oontent  ourselfes  with  methods  which  satisfy  tho  following 
conditions ; 

1. )  the  first  node  (a^»o  )  has  multiplicity  m  (m£2)  ; 

2. )  all  further  nodes  have  multiplicity  j  2  | 

thus  (cf.  (1.18)  ): 

(3.1c)  (is2,...,n  >  k=3»4>*«»»ifl) 

(3. Id)  b^j^=o  ( j“2 , • • * , i-1  ;  k=3» « « * »m) 

3. )  m  and  s  satisfy  tho  condition 

(3.2)  m+s  $2m+2  s4'm+2 


Conditions  for  the  Coefficients 


Satisfaction  of  C(m) 


0=1  0!  j=1 


) 


ij  5 


^  (krJ  i-.li.q.-l)¥~b  .^)a^-q+ 
1  >2  ij  J 
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1o9 


+  ib(kU 

k  il  k 


( i—  1 1  .  •  •  |  n  I  k*  1 1 1 1 1 1  n ) 


For  i=1  and  k*1,...,m  this  is  satisfiod  because  of  a^>0| 
by  putting 


(3.3a) 

(3.3b) 


( k=  1 1  •  •  •  y  n) 


(i-3 . .  ,  k-1 . m) 


C(rn)  i3  satisfied  completely  and  the  coefficients  a^t,,,tan  , 
b| (is3f . . .  |  j® 2 , . . . ,  i-1  |  1  f  2  )  arc  still  free • 


Satisfaction  of  B(m)  : 

B(m)<f=r>  ^  (k-1i-1JC-l)f"c^c)aJ"CJ  (Vn) 


on 

k 


fn 


Jl(k-1i-1»C-l)57c^^n^"°  +(k-1 ) !  e *  £  (k*1,...,m)  | 

C=1  1*2 


by  putting 


(3.4)  *  j.’tTTj!  ^  (k”1  >*1  »°*1 )  J^°^g^ak“g(k-1>...  >m) 


,00 


B(m)  is  satisfied  and  the  coefficients  '  (i*2,...,n  ;  k*1,2  ) 
are  still  fret. 


Using  (3.2)  the  elementary  weights  can  be  distingished  as  follows: 

1 .  )  l  =  C  fy  *  *  with  m  (i*1 , . . .  ,t)  : 

because  of  C(m)  Lemma  1  (  from  seotion  V.2,  )  is  applicable 

I*-  1  n 

and  all  these  weights  can  be  reduced  to  J  (by  (2.12)  ). 

2. )  0  =  with  r1^  m  f  ri^  111  (i=2, ...  ft)  : 

all  these  can  by  Lemma  2  (  cf.  (2.16)  )  be  reduced  to 

”,■’3  . 

The  condition  (3.2) 
m+s  occur  in  the 


guarantees  that  all  elementary  weights  of  order 
above  cases  (  cf.  plol  ,  2.)  ) 


11o 


V.3. 


We  thus  oan  restriot  us  to  the  elementary  weights  of  the  forms 

[Z1**1]  and  with  r  i  tn+s  ,  r^  >  m  ,  k  . 

Because  B(m)  and  the  Corollary  of  Lemma  1  (  in  V.2.)  by  (3»3&|b)  and 
(3*4)  already  all  conditions  for  the  orders  £m  are  satisfied.  We  thus 
oan  restriot  us  to  the  rest  and  determine  the  still  free  coefficients 

a2,***,an  }  0i^  *  bij^  (i“2#...,n  ;  j»2,...,i-1  \  k-1,2  )  to  satisfy 

the  conditions  for  the  orders  m+1,...,ra+s  . 

V/e  thus  obtain  the  following  conditions  for  the  coefficients: 


-rder  m+1  : 
r/n-f  1 

L*  J  *  ST  1 


order  m+2 
/m+1 


1 

m+2 


[p  ] 

nyjj  •  ffTTTffgy 


order  m+3  : 


li/B+1J].  (b+2j}.+5)  lOT/U  -  T5TTO+F7 

Generally  the  additional  conditions  for  order  m+s  (s^  m+2  jof.  (3,2  )) 
are  as  follows  : 

i,...lwZ  JZ  -'•••Z  ->  “  (m+k^  +1 )  (m+k^+kg+2  7777  (m+k^  +  .  *  T+iTTrJ 

with 

k.+k0+, ,,+k  +t»8  whore  k.  £  o  and  T-1t»»»f3  • 

1  t  T  1 

.o  can  be  seen  oasly  (induction  !  )  for  order  m+s  there  exist  2s 
oonditions  of  this  form. 

’Jp  to  m+5  these  are  listed  in  /2 5/ » P -45  • 


V.3. 
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The  method  of  Fehlberg  as  special  oase 


a)  Fehlberg1 s  method  which, is  known  for  its  acouraoy  and  its  low 
expenditive  of  work,  works  as  follows  (cf,  / 1 2/ »/l 3 /  ): 


diff.  equation:  y'  *  f(x,y)  .  y(xQ)  »  yQ 


approximate  solution  of  order  m:  y(x)=  y]0  “  y  h*  YL 

k»o  k«o 


diff.  equation  for  the  approx,  solution:  A 

$’(x)  -  f(x) 


khk"1Y,. 


V/o  put:  z(x)«  y(x)-y(x)  ,  where  y(x)  is  the  exact  solution  and 
obtain  for  z(x)  the  following  diff.  equation  : 

(3.9a)  z'  (x)=y'  (x)-y’  (x)=*  f  (x,y(x)+z(x))  -f(x)  «:  f(x,z(x)) 

with  z(x  )»z'(x  )».,.«z^(x  )«o  . 

'  o  '  o'  o 

This  we  now  solve  by  a  Runge-Kutta-process 
n 

(3.9b)  z(x)  =  h  ^ciki 

i-1 

(3.9c)  -  f^+a^^+h  ^b^k.)  (i«2,...,n)  , 

for  which,  because  of  the  above  made  transformation,  considerably 
higher  orders  are  possible: 
improved  solution  of  order  m+s: 

/  \ 


n 


(3.9d)  y.(x)=  y(x)  +  z(x)  =  y(x)  +  h  21  o.k. 

1  IT2  1  1 


Remarks : 

1.  In  / 1 3/  Fehlberg  gives  a  method  of  order  m+4  using  six  nodes  and 
with  an  estimation  of  the  leading  error  torm. 

2.  In  /5l/,p.1o1  Y/anner  has  extended  all  the  theory  of  Butcher  to 
Fehiberg-processes . 

b)  Rungc-Kutta-processes  with  one  m-fold  node  : 


Here  we  are  interesting  in  process  with  one  m-fold  node  and  a  f> 
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sf'C'n-l,  single  nodes  . 

Let  this  fact  be  expressed  by  the  following  diagramm  : 


m  o 
o 


2  o 

1  4-o~o— o- 


1  2  3 


n 


The  nodes  of  the  method  are  piotured  on  the  abscissa  and  the  ordinates 
show  their  multiplicity. Prom  this  diagramm  it  is  now  olefr  why  in  (3.9V 
and  (3«9o)  the  indices  i,j  start  with  2. 

Both  methods  are  now  shown  to  be  identical: 
no  have  from  (1.15*16)  *• 

(Mo.)  y2W-  yQ+h  |^(1)+h2o1(s)g1(2)+...+h'"c<n)g1("> 


no  have 

from 

(3.1oa) 

y2(x 

with 

(3. lob) 

g(k) 

61 

(3. loo) 

m 


(i-2,... fn) 


and  it  holds  the  following  theorem: 


whoorQT :  The  method  determined  by  (3*1oafb,c)  is  identical  with 
Fchlberg^  process  (3.9a-d)  if  we-  put 

(3 . 1  od )  cP^-c.^  (i-2,...tn) 

loe)  "fte'lT) ;  "^0i‘i)ai~1j  (k-1 , . . . ,ra) 

bii 

ij  ij 

(3. log)  (i*2 * • . •  *n  ;  k*1»...»m)  , 


(3. 


O.lof)  Ta^X  b.  .  (i*2 , . . .  ,n  j  j-2,.,.,1-1) 


Proof :  We  have  to  show  that  y^(x)  «=  yg(x): 
inserting  (3.1  >b,c,e)  into  (3»1oa)  wc  obtain 


j 
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y2(x)  -  *\-£khk  ^o^af'vh  . 


?(x)+h  j- 


9(x)+h  (xQ+a^jj 


(3.11a) 


y(x)+h 

with 


(3.11b)  kx  -  g|1^-f(xo+aih)  (i«.?,...,n) 

Next  we  insert  (3.1ob,f,g)  into  (3. loo): 

fi(l)»  f(x  +a  h  v  +hb  ) a ( 1 )  +  ^  ) )  ■ 

Ri  "  f^xo+aift,yo  hbi1  g1  +,,,^Ti1  g1  +h  2-®iJ  ej  ' 

f  (xo+aih,yo+aihY1+*  •  ,  +  (aih)  Ym-^kh  ^bij^nj  Yk+h^bi,Vg,1  * ' ; 
f  (xo+aih*y(xo+aih)+h^  b}^(gj^-^k(a.jh)k“1Ylc))  ^5*-1b^ 

(3.Ho)  f(xo+aih,^(xo+aih)+hV^b^k^)  (i«2,...,n) 


Prom  (3«Hb,c)  vc  have 


(3.12a)  ki  =  f (xo+ajh,y(xo+aih)+h^b^kJ)  -f^+ajh)  (i-2,...,n)  ; 

J=2 

On  the  other  side  it  follows  from  (3*9a,c)  • 

(3.12b)  kA  «  f (xo+aih,^(xo+aih)+h^bijkj)-f(xQ+aih)  (i-2,.,.,n)  . 

Next  compare  the  conditions  which  are  to  bo  satisfied  by  ^  and 

b^j'  (i-2,...,n  j  j»2,...,i-1)  with  those  for  the  method  of 

Fehlborg  (e.g.  in  the  form  given  by  Wanner  /5l/»P»1o3  )•  We  therebj 
confirm  easily  that  the  equations  coincide;  thus  we  may  put 

(D. 


bi^"bij 


and  *0^  (  ia2 y • • • pn  $  j*2 f  •  •  •  y  i*1 )  • 


Hence  from  (5.l2afb)  we  have  k^-k^  (i**2,...,n)  and  (3.1  la) ,  (3«9d)  shows 
that  y«(x)-y2(x). 

The  coefficients  (3.1oo,6)  aro  of  course  determined  by  (3*3)#(3«4)  • 
Done. 


Some  explicit  methods  of  orders  (m+2 ) , . . . , (ra+5) 


Hero  we  list  some  explicit  formulas  of  different  orders.  Their  deri¬ 
vation  from  the  above  conditions  is  given  in  full  detail 
in  the  thcsic  / 2 ,  pp  51-61  . 


1 . )  Formula  of  order  m+2  : 

n 


diagrsn 


2 

1 


<1  o 


c— -e 

1  2 


n*»2 

m  >2 

order  m+2 


Coefficients : 
ar°  ’  V  m+2 


■  0).i  .0). 


=1  ,  c; 


3(2)=  _ _ 1 _  . 

\"‘(r.i+1 )  (m+2  ) 


ffc)  1  !  1  /.  (2)  k-2\  0 

!1  -  1F-T;T  ik  ■  (k-l)c2  *2  )  (kl2 . . 


u 


( k”  1  > .  •  • » 'l ) 


No  coefficient  oan  be  choosod  freely. 


2 , )  Formula  of  order  n+5 : 

i 

m 


diagram 


2  <•>  ?  o 

1  o-o—o — - 

1  2  5 


n-2 
m \2 

order  m+3 


oocffiy ionto  i 


m 


v°  >  V  ST}  ’  V  1 

0^-1  ,  I 


,(2). 


,00. 


3(m+1 )(m+2)a 


2 _  (2)  1 

m  *  3  3{m+l ) (m+2 )  * 


|kliy;'  -  (k-1)(c^ag“2+c^  )j  (k-2 , . . .  ,m)  { 


^2  (k“ 1 » • • • »m)  j 


V,  v  '  /_n 

52  ’  52  ' 

v(k)  k(k-1)  (2)  k-2 

b3l  1  2  b32  a2 

( k~  1 , .  #  •  f  ra) 

• 

3.)  Formula  of  order  m+4 

• 

• 

m 

i 

<D 

diagram:  2 

1 

i 

4 

(j  0  0  0 

t 

#1  t\  r\ 

n-4 

o*2 

order  m+4 

1  2  3  A 

Coefficient  a  : 

V°  »  V  m+6  »  V 

n+2  „ 

“~r  •  a  .  ■  1 

m+4  <+ 

} 

.0) 


,0).  .(i)_  „0). 


or7-i . °rao* 'bo 


(m+6) 


3  "w4 

2 


(2).  _ 


,(*) _ _  c  „ _ 

12a®  (m+1 )  (m+2 )  (m+3  )2  ’  3  4a®(m+l)(m+3)2 


,(2)  m 

4  “  Ffm+1  )  (m+2 )  (m+3 ) 
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,(l)_  (2)_  _ 2 _ 

32  ’  32  c (2  ^  (m+1 ) (m+2 ) (n+3 ) (m+6) 


b(2)=  _ _ 1 _ 

Tp^a^^in+I )  (n+2  )  (rn+3 ) 


-  b<2>- 

2  *  b43 


cjj2  (m+1 )  (m+2 )  (m+3  )2 


,  (k)  .k  kfk-OTT4,  (2 )  k-2  /.  ,  .  ,  „  \ 

bi1  "  "i  -  -V-iAb ij  “j  f1-5-4-  •  ** . .  > 


4 . )  Formula  of  order  m+5  : 


diagram : 


0  O  0  > 


1  2  3  ••!  5 


n«*  5 

m  »3 

ordor  m+5  • 


Coefficients : 


ar°  »  v  i+j  *  "3“  m+3  »  v  i+5 » nr 1  * 
oS1).,  ,  =  0>.o0>.o  , 

0(2) _ Jsl _ 0(2). _ i. _ 

loa^*2  (m+1 )  (n+2 )  (m+3 )  (n+4 )  3  6a™+2 (n+1 ) (m+2 ) (m+4) 


n  _  £+2  =1  . 

n,r*  mxi;  »  ae;  1  i 


A’  6ar?+2  (m+1 )  (m+3 )  (m+4 ) 


(2)_  3m  +1 5m+1o 

’  0  5  1 5 (m+1 ) (m+2 ) (m+3) (m+4)  } 


°i(k)=  *  (*-0(^2)^...^2^-2)j  (k. 

b2l^=  a2  (kss1  >•••♦*) 


2,..., in)  } 


V1*-  b{2)=  o  • 

\32  '  b52  °  * 

’0(^)_  b^^-o  b^2^=  t' -v - — i  — 

42  '  *43  ’  V  7T2j.nk  iw^. 


3c;  'n£  fm+1 )  (m+2 )  (n+3 )  (m+4 ) 


b(2)-  o  - 

'•3  7 


b(0.  b0)  b0)_  0 

ucp  -  cr*  nr./  “  0  » 


"52  53  ”  54  “  ’ 

v  (2 _ J. _  v(2). 

Ccp^a^  (m+1 '  (n+2  )  (n+3 )  (rr+4)  ^ 


m+3  _ 

2c^2  m(m+1 )  (m+2 )  (i..+ 


v.3. 
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b 


00. 

il 


1 


2c 


T27-m 


a"(ra+1)(ra+3)(m+4) 


i-1 


Qk  _  k^  f  b(2)aK 

i  2  ^5°tJ  aj 


k-2 


(k-l,.,.,ra  ;  i«3,4,5  ) 


Nunerioal  Examples 


The  practical  evaluation  of  these  formulas  is  only  valuable,  when 
it  is  combined'  with  the  use  of  the  recursion  formulas  which  are  descri¬ 
bed  in  Chapter  II.  With  these  the  methods  again  can  fully  be  made  auto¬ 
matic  by  using  the  3ame  subroutines.  It  nay  finally  be  noted  that  in 

2 

many  cases  the  calculation  of  Df,  D  f,...  often  requires  much  less  work 
than  the  calculation  of  f  itself  (ospecial  if  there  occur  a  lot  of 
olomontary  functions  like  exp,  log,  sin,  cos,  ...  ). 


In  the  following  the  methods  arc  tested  at  some  differential  equations 
with  known  solution.  They  are  further  compared  with  the  method  of  Fohl- 
borg  and  with  the  powor-sorios  method.  All  computations  are  with  order 
10  and  have  been  carried  out  with  single  precisions  (9D)  on  the  Zuse  Z23 
computer. 

-x2  \ 

Example  1  :  y'  =  2x(e  -y)  y(o)*=1  » 

?  -*2 

solution:  y(x)=  (1+x  )o  , 

Example.  2  :  y'=  -|y2x“^2  y(l)=1  $ 

solution:  y(x)=  Nx  • 

Example  3  :  y '=  1-e"^(sin  x  -  cos  x)  y(o)=o  j 
solution  :  y(x)=  log(sin  x+ex)  . 


example  4 


y'  =  cos  x  .(y+sin  x)  y(o)=>1 
solution:  y(x)»  2osin  x-sinx  -1  . 


Example  5  :  y'  =  (x4+y4)/xy5  y(1  )=  1  , 

solution:  y(x)*  x. (1+4. log  x)1/^ 
Example  6  :  y*  =  2(xy^2-y)  y(>)~  o,25  » 

solution:  y(x)=  (e  +x+l)“ 

Example  7  :  y'=  (xy  +y)/(x.log  x)  y(o)=  o,5 

solution:  y(x)=  log  x/(e+2-x) 


In  the  following  table  the  errors  of  the  different  methods  with  those 
examples  and  with  the  given  stop  sizes  are  listed. 


Example 

powor 

h  jjsories  j  I 

“-L  JILL 

Fc 1 

1. 

o,5 

1 ,6. lo”6 

1,6.1o’7 

4,3. lo’8 

1 ,5.1o”6 

8 ,2 . lo’8 

n 

. 

2. 

o,8 

4,o.1o’^ 

5,4. lo’6 

8,1,1 o"7 

4,4. lo”8 

7,6. In’8 

6  , 4 .  * 

3. 

o,25 

5,2. lo’6 

1,o.1o’7 

4,4. lo”8  | 

8,5.1o”10 

1,o. lo’8 

5  ,  •  1 

4. 

1 

i 

o 

• 

ON 

3,4. lo’6 

2,8.1 o’6 

3,2. lo’6 

i  i  "6 

1 ,0. lo 

1,3.1 

5. 

o,2 

1,4. lo"3 

l,5.1o*5 

4, 2.1  o’6 

3,7. lo”7 

4,8. lo’7 

1,1.1 

C . 

0,4 

9,2. lo’6 

5,2.1 o’7 

1 ,9.1o"7 

2,4. lo’8 

4,3. lo’8 

1.0,1 

7. 

o,7 

1,1. lo"5 

1,0. lo’6 

1  ,o.1o’6 

7,7.1o’8 

4,2.1o’7 

1,-1 

h=  step  size 

I  :  formula  of  order  in+2  j 

II  :  formula  of  order  m+3  \ 

III  :  formula  of  order  m+4 

IV  :  formula  of  order  m+5  . 

As  can  bo  seen,  the  rosults  of  formulas  III  and  IV  on  the  average  have 
the  same  accuraoy  than  the  method  of  Fchlberg.  It  can  further  bo  seen, 
that  with  equal  order  the  methods  with  more  nodes  are  very  much  better. 
The  results  of  III,  IV  and  Fohlberg  arc  mostly  2-3  digits  better  than 
with  the  power  series  method  of  the  same  order.  In  addition  note  that 


the  necessary  work  for  the  power  series  method  mostly  is  higher  than  with 
tho  other  formulas. 

The  following  questions  are  still  open: 

1.  "  optimal"  methods:  tho  coefficients  in  general  are  not  uniquely 
determined  by  the  conditions.  Some  of  them  have  been  fixed  ar¬ 
bitrary,  mostly  to  reach  simple  results. 

How  are  they  to  be  fixed  to  give  methods  with  minimal  error? 

2.  How  arc  effective  error  estimates  possible  ? 

3.  How  is  the  stability  of  the  methods  ? 
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Chapter  VI 


On  Step-size  Control 


by  G.  Wanner 


This  chapter  dcal3  with  the  problem  of  choosing  the  step  sizes 
in  the  numerical  integration  of  ordinary  differential  equation  . 
using  one-step  methods.  First  the  frequently  used  formulas  arc- 
discusse-d  which  try  to  keep  the  local  error  constant.  Then 
expressions  for  an  ’’optimal”  step-size  control  are  developed 
which  take  into  account  the  propagation  of  the  local  errors 
to  the  final  result.  Numerical  results  are  given  and  compared 
with  those  of  the  step-size  control  of  Morrison. 


VI.  2. 
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VI. 1  dtep-size  Control 

A  system  of  n  ordinary  differential  equations 
y*  =  f(x,y) 

is  given  and  the  solution  y(x)  with  initial  values  xQ,y0  is 
wanted  at  some  point  x^.  Using  some  one-step  method,  the  inte¬ 
gration  proceeds  on  the  steps  x0<xi<x2<# • *<xtf  with  the  step- 
sizes  h1=x1-x0,  h2=x2~x1, . . .  jh^sxjj-Xjj^j.  For  a  step  size  control; 
the  method  has  to  be  equipped  with  some  error  estimation,  i.e., 
to  each  initial  point  x^.y^  and  step-3ize  hk+1  it  gives  a  approx 
mat  ion  to  the  solutions  and  approximate  error  estir.atio 

R^.  The  usual  procedure  now  is  trying  to  keep  these  local  errors 
equal  to  some  0iven  numbers  y^,  the  wanted  errors.  These  might  o 
lo-5,lo~10,lo~20  and  depend  on  the  wanted  accuracy.  Thus  by 
putting  !R  : 

(1.1)  n  =  max----  - 

i  Yi 

one  tries  to  keep  nHl.  A  possible  procedure  is  noi>;  the  following 
The  first  step  is  calculrr  .Uvi  a  guessed  ,rtep-size  h^.  Then  n 
can  bo  evaluated  by  (1.1).  Of  course,  n  will  not  be  equal  1.  If 
p  is  the  order  of  the  method,  a  much  more  better  step -size  would 
have  been 

(1.2)  h  =  h1ei^Vl7n. 

But  if  n  ia  not  /ery  much  greater  than  1,  say  n<n2  with 
n2s  1.5  or  10,  then  we  use  h  for  the  next  step 

ho  “h  , 

4ST-. 

otherwise  we  re  peat  the  first  step  with  the  step  size  h*h^. 

The  same  procedure  is  then  also  used  in  the  following  steps 


VI . 2 .  Damping 

Occasionally,  especially  in  regions  where  R^  changes  sign,  n 
may  be  very  small  or  even  zero.  In  such  cases,  formula  (1.2) 
would  lead  to  an  excessive  increase  of  the  step  size.  For  this 
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reason,  one  chooses  a  number  n1<l  (say  sl/10,  1/100)  and  if 
one  replaces  (1.2)  by 

(2.D  e  =  hj  m  a  -  af 

Thus  the  step  size  can  increase  at  most  by  the  factor 


This  stabilizes  the  step  size  control  and  guards  against  over 
flow.  Formula  (2.1)  is  obtained  by  replacing  the  hyperbola 
P+1  / — ' 

vi/n  by  its  tangent  of  the  point  n1* 


VI. 3.  ilorrison’s  Control 


Consider  the  example 

(3.1)  y*  =  y2,  y(0)=l,  y(0. 999999)=? 

The  solution  is  y=l/(l-x);,  y (0.999999) =10^.  The  solution  for  a 
general  initial  value  y0  is  y(x,y0)  =  l/(l/y0~x).  The  derivative 
of  this  solution  with  respect  to  that  initial  value  is 


(3.2)  H(x)  =  =  l/(l/y  -x)2y2  =  y2/y2. 

3y0 


Thus,  if  the  initial  value  y  =1  is  changed,  say,  by  10  *#,  then 

0  -3 

the  solution  at  the  point  0.999999  changes  by  10  J  since 
12 

H(0.999999)=10  .  If  we  compute  this  example  by  a  stepuise 

-19 

numerical  integration  with  local  accuracy  10  %  the  final 

result  will  not  be  better  than  10*  \  Of  course  here  it  is 
unwise  to  compute  also  the  last  ctep3  with  this  same  accuracy. 
The  last  steps  need  not  to  bo  calculated  with  the  same 
accuracy  than  the  first.  The  idea  lies  at  hand  to  multiply  the 
chosen  error  sizes  y^  by  the  connection  matrix  H(x)  along  the 
solution.  This  means  to  replace  (1)  by 
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(3.3)  n  =  mgx|  • 

This  is  the  step  size  control,  which  Morrison  /37 /  has  proved 
to  be  nearly  ’’optimal'*  for  the  case  when  n=l  (one  equation 
only)  and  when  the  errors  all  have  the  same  sign. 


VI. 4.  Another  Possibility 


There  is  still  another  possibility  for  a  step  size  control 
which  shall  be  derived  now: 

Assume  the  differential  equation  to  be  integrated  from  xQ  to 
xN  using  H  steps  hjsx^-x  ,  b2=x2”xl*  *  *  *  ‘  Thc  local  terror  of 
the  j-th  step  wc  denote  by  and  its  propagation  to  the  final 
result  is 


(4.1) 


e!N^  =  H(xN)H~1(x^)eJ.  . 


We  again  assume  that  n=l  anu  that  all  errors  have  constant  sign 
although  the  results  may  bo  interpreted  for  the  other  cases 
as  well. 

Neglecting  rounding  errors  we  nay  assume  c.  =  <f>,h?+  ,  thus 

v  v  J 

e!N)=  Xjh?+1  where  x^sH(xn)H“1(Xj  )^. . 

Assuming  that  does  not  depend  on  hj,..,h.  (what  however 
actually  is  the  case),  we  solve  the  easy  minimum  problem 


N 

J=1  13 


N 


■ao  =  T~  xjhj 
J=i 

under  the  condition  that 

•1 


P+1 


=  min! 


\  h.  =  x..-x_  =  C  . 
1=1  d  l*  0  0 


The  method  of  Lagrange  multipliers  gives  h.  =  C1/  W-  ■ 

J  ^  J 

Thus  the  local  error  e,  should  be 

J 
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vvT1" 


k  h 


*.iC2 


3 


icnce,  the  sten-size3  arc  chosen  “optimal11 .  if 


(4.2) 


ten- 

H(x..)lf*(x,)e.h71  =  k, 

N  J  J  J  /  M \ 

i.c.;  if* the  contribution  ej  '  of  each  step  to  the  final  result 

is  proportional  to  the  step-size . 


In  the  case  n=l  H(xN)  is  only  a  constant  number  and  need  not  be 
knov/n.  Hence,  in  the  course  of  computation,  it  is  only  necessary 
to  keep 

(4.3)  H"1(x.)ejhT1«Y« 

This  result  differs  from  that  of  Morrison  by  the  denominator  hj . 
An  error  estimation  for  the  total  truncation  error  is  now  ob¬ 
tained  as  follows: 

The  final  error  ej1^  (4.1),  which  results  from  e^  is  because  of 

(4.3)  equal  to  *hjH(xN)y.  These  numbers  are  summed  up  easily  to 
give 

(4.4)  E  =  (xN-x0)H(xNb 

as  estimation  of  the  total  truncation  error  of  the  final  result. 


Formula  (4.2)  may  also  be  interpreted  for  systems  of  differential 
-quations.  But  then  the  knowledge  of  the  final  connectionmatrix 
H(xn)  is  necessary.  On  the  other  hand,  the  use  of  (4,3)secms  only 
adequate,  of  all  eigenvalues  of  II ( x^ )  have  approximately  the  same 
3ize.  This  is,  because  (4.3)  nappo  the  error  not  to  the  endpoint 
xj(JJ  but  to  the  initial  point  xQ  and  the  box  maxiRj|3y  may  chanrv'-' 
shape  considerably  11 
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Finally  v;o  mention  the  paper  Greonspan-Hafner-Ribaric  /20/. 
There  Morrison’s  control  is  compared  with  several  others 
(such  as  the  ’’natural”  step  size  of  Collatz  /7/,p.89).  For 
differential  equations  with  constant  coefficients  (y'scy) 
Morrison’s  control  as  well  as  the  above  ‘'optimal”  control 
become  a  control  with  constant,  step  size-  (for  the  methods 
considered  in  this  report). 


VI. 5.  Numerical  Examples 


Uainc  the  Lie- series  method  of  Chapter  III,  the  different  con¬ 
trols  of  the  sections  VI. 1,  VI. 3  and  VI. 4  have  been  compared 
at  several  examples.  The  results  are  listed  below: 


Example  1;  y’ 

=  y 

2,  y(o)  =  l 

solution 

y(x)  =  1/(1 

-x) 

step  • 

actual 

error 

time  P'- : 

size 

• 

error  of 

estimation  1 

3tC-p 

control  m Is 

„  1 

V 

Y  '  X 

y(x) 

;  for  y (x) 

steps 

(m  sec' 

normal  10 1 3 

2 

10“ 

^•0.9 

3.01*10~ll> 

:  — 

23 

44 

!  1 
i 

!  1 

| 

j 

i 

0.99 

3.05*10" 13 

i  .....  I 

' 

50 

' 

I 

! 

0.999999 

2.89*1C~5 

l 

j 

210 

15  :3 

2 

10" 

19 i 0.9 

1.65*10'  17 

J 

j 

15 

60 

; 

i 

0.99 

1.67* 10“15 

1 

I 

32 

i 

; 

0.999999 

1.73*1  O'7 

123 

optimal  10  3 

r  • 

2 

10“ 

-i '/  ’! 

1  0.9 

'6.8*  10~16 

•  9.0*  10*  | 

25 

56 

1 

0.99 

:  7. 5*10“^ 

-1 4 

j  9.9*10  iq, 

50 

0.999.999 

7.6*10'  b 

,10.0* 10“& 

148 

19  3 

2 

10 

^^0.9 

6.4* lO"12 

i  9.0* 10  18 

16 

71 

j 

> 

;  0. 99 

.  —  1 6 
7.0*10 

9.9*  10“'16 

31 

;  i 

i 

0.999999 

7.2* 10“ 8 

10.0*10  8 

92 

Un  to  0.99  onl y.  the  optimal  step  size  control  fives  no 
noticeable  increase  of  effectiveness. 
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t  » 

Example  2;  y^y^  y^y^,  Y1(0)=0,  y2(0)=l,  solutions:  y^^  =  sinh 

y2  =  cosh 

Results  for  x=10: 


step 

size 

control 

m 

i 

B 

l 

k 

1 

1  Y 

j  actual 
| error  of 
h  j  y±(10) 

error 
estimat. 
for  y. 

steps 

normal 

13 

5 

it 

0.83  ..0.53  jl.21*10"i6 

15 

optimal 

13  i 

5 

3 

!io-20 

0.805  ..0.811  |3.75*10“16 

4.4*10“16 

13 

optimal 

13 

5 

L_ 

3 

! 

10~21 

0.7203. .0.720514.2  *10-17 

4.9*10'*1V! 

14 

As  expected  for  linear  systems  with  constant  coefficients,  the 
optimal  step  size  remained  constant.  The  estimate  for  the 
total  propagated  error  is  satisfactory. 


Example  3:  y ' =-xy3 


y(-U  *  y0. 


solution  y(x) 


yo(l+(x2~l)y2) 


a)  yQ  =  0.999999995,  thus  y(0)  *  10000 

y(l)  =  0.999999995. 

Results  for  m=15,  5=3,  k=2: 


step 

size 

control 

Y 

h 

actual 
i  error  of 
x|y(x) 

error 
estimat. 
for  y(x) 

steps 

normal  * 

0.2 

...  0.00003 

0;1.9*10fl 

— 

28 

j 

i !  i .  s^io"*11 

— — 

67 

optimal 

!  -12 
10 

0.2 

...  0.00007 

0  4.9'ID"1 

Q.Ono"1 

21 

< 

1 14.1*10"15 j 

20.0*10'13 

J 

55 

For  this  example •  constant  stop  size  is  not  advisable. 
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b)  yQ  =  0.99999999995,  thus  y(0)  *  100000 

y ( 1 )  =  0.99999999995. 

Results  for  111=15,  s=3>  k=2: 


step 

size 

control 

Y 

h  i 

t 

X 

jactual 
icrror  of  i 
y(x) 

error 
estimat . 
for  y(x) 

steps 

normal 

10"12 

0.29. .0.0000023 

0 

1.6  10+5  j 

— 

46 

1 

1.7  10-13 

j 

107 

1 

optimal 

io"w 

0.24. .0.0000060 

1 

0 

5.1  10“° 

9.6  10*° 

32 

1  .. 

|l 

1 _ 

5.4  10“15< 
- - - 1 

20.0  10* 15 

i  80 

* 
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Chapter  VII 


\  * 

Calculation  of  Switch-on  Transients  at  the 
telegraphic  Equation 
hy  Means  of  Generalized  LIE  -  Series 

*  __  _y** 

hy  R.  Saoly 


Abstract 


This  ohapter  deals  with  the  switoh-on  transients  occurring  in 
the  telegraphic  equation, i.  e.,  an  initial  and  boundary  value  pro¬ 
blem  of  a  hyperbolic  partial  differential  equation,  by  means  of 
generalized  Lie  series.  \!q  shall  assume  that  an  ordinary  alter¬ 
nating  voltage  U(o,t)  -  A  coswt  +  3  sinojt  is  applied  across  the 
input  terminals  of  a  telegraphic  line  (electric  twin  lino)  of 
length  a.  V/e  confine  our  investigations  to  two  limiting  cases, 
namely,  that  the  line  is  either  shorted  or  open  at  the  other  end. 

The  first  part  of  the  paper  gives  a  formal  solution  using 
power  series.  The  solution  is  represented  by  moans  of  Lie  scries 
with  n  generalized  Lie  -  Operator.  Next  the  switch-on  transients 
is  treated  with  shorted  wires  and  given  initial  and  boundary  con¬ 
ditions.  Two  numerical  examples  shall  illustrate  this  switch-on 
transients  problem.  Finally  the  computation  of  the  solution  U(x,t) 
and  j(x,t)  for  the  initial  and  boundary  value  problem  with  open 
wires  is  given. 

My  thank  go  to  Prof.  W.Groebncr  for  his  assistanoojl  wish  ack¬ 
nowledge  tho  discussions  with  H.P.oitbergcr,  G. Wanner  and  K.H.  Kast- 
lungor. 
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VII .  1  Introduction 


1.  The  Telegraphic  Equation 


In  the  present  paper  wg  shall  calculate  the  switch-on  transients 
occurring  in  the  telegraphic  equation,  i.  e.f  an  initial  and  boun¬ 
dary  value  problem  of  a  hyperbolic  partial  differential  equation, 
by  means  of  generalized  Lie-series.  We  shall  assume  that  an  ordinary 
alternating  voltage  U(o,t)  =  A  cos  wt  +  B  sin  wt  is  applied  across 
the*  input  terminals  cf  a  telegraphic  line  (electric  twin  line)  of 
length  a.  We  confine  our  investigations  to  two  limiting  cases,  namely, 
that  the  line  is  either  shorted  or  open  at  the  -'ther  end. 

Lot  a,,  a2  bo  tho  input  and  o, ,  ^  thc  output  termlnals  of  tho 

line.  V Jo  take  one  axis  of  coordinates  along  the  line  and  denote  the 
distance  from  the  input  terminals  by  x.  The  length  of  the  line  is  a. 

At  time  t,  the  current  j(x,t)  flows  in  the  wire  at  the  point  x,  the 
voltage  between  thc  two  wires  is  TT(x,t), 


a. 


*j(x,t) 


\V 


x^=o  x  x=  a 

Fig.  1.  Schematic  diagram  of  a  telegraph  line 


Of  these  parallel  wires  wo  consider  a  very  small  line  element 
I, II, III, IV  of  length  dx  (infinitesimal  four-pole),  w©  assume  the 
lino  constants,  referred  to  unit  length, to  be  independent  of  space 

j 

and  time' coordinates  and  denote  them  by  the  following  symbols: 


r  .  resistance 

1  .  inductive  reactance 

q  .  conductance  (leaking  insulation) 

c  .  capacitive  reactance 
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r  dx 


1  dx 


IIP 


aJr  1"  1  ?7?7WV- 

_ S 

IT  W  ■  l  |  L  [  ' 

\  - 

o  dx 

q  dx 

1  n 

' _  _ _J 

L 

dx 


J  4  Jxdx 
U  +  Uxdx 

IV 


Pig.  2.  Infinitesimal  element  of  a  telegraphic  line 


The  two  basic  equations  of  the  telegraphic  equation  follow 
from  the  laws  of  electromagnetic  theory: 

(1.1)  JxU,t)  =  -QU(x,t)'-  cUt(x,t) 

(1.2)  Ux(x,t)  -  -rJ(x,t)  -  lJt(x,t) 

To  eliminate  current  in  (1.2)  wo  differentiate  the  first  equatio 
with  respeot  to  t  and  the  second  with  respect  to  x.  Putting 

tq  *  a 
rc  +  q1  ■  (3 
lc  =»  y 

re  obtain  the  telegraphic  equation  for  the  voltage  U(x,t): 

(1>3)  Uxx(x,t)  ”  aU(x»t)  +  +  yUtt^x,t^ 

An  analogous  equation  in  J(x,t)  oan  be  found  by  differentiating 
(1.1)  with  respect  to  x  and  eliminating  Ux(x,t)  by  moans  of  (1.2): 

(1.4)  Jxx(x,t)  =  ccJ(x»'fc)  +  PJt^x»t)  +  yJtt^x,t^ 

In  their  physioal  meaning,  the  constants  a,  (3  and  y  are  positive 
’^thematical  treatment  requires  only  that  y>>o,  beoauso  this  makes 
the  equation  hyperbolic.  With  a  =  (3  =  o,  the  telegraphic  equation 
becomes  the  ordinary  wave  equation. 


Vli'.l . 
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2.  Formal  Solution  of  the  Telcgraphio  Equation 


Equation  (1.5)  describee  all  electric  phenomena  in  the  two 
parallel  wires.  A  problem  frequently  arising  is:  what  happens  when 
the  system  is  switched  on?  At  given  time  (t=o)  the  system,  whoso 
electric  condition  at  that  moment  i3  known,  experiences  some  kind 
of  external  influence.  We  wish  to  calculate  the  changes  produoed  by 
this  influence.  Wc  assume  this  influence  to  be  a  suddenly  applied 
voltage.  We  assume  further  that,  at  the  point  x»o,  the  voltago 
U(o,t)  -  f(t)  is  a  given  function  of  time.  From  the  basio  equation 
(1.2)  we  find 

Ux(o,t)  -  -rJ(o ,t)  -  Ut(o,t)  =  -  rg(t)  -  lgt(t)  «  h(t ) 

9 

Now  we  shall  solve  the  telegraphic  equation  with  the  power 
series  expansion  (cf.  /2j/  ,  p.112) 


(1.5)  U(x,t)  -2H 

V=*  0 

where  the  functions  ( t )  arc  yet  to  be  determined.  The  functions 

iQ(t)  and  i'^(t)  can  be  found  from  the  initial  conditions: 

(1.6)  u(o,t)  »  f(t)  =.  ¥  (t)  and 

(1.7)  Ux(o,t)  -  h(t)  -  ^(t) 

be 

The  remaining  functions  V  (t)  mayVralculated  by  means  of  a  re¬ 
cursion  formula  which  can  bo  obtained  by  comparing  coefficients  of  xv 
after  the  power  series  expansion  has  been  inserted  in  (1.5). 


(1.8) 


(v  +1  )  (v+2 ) 


All  functions  'Fy(t)  can  be  calculated  from  this  formula,  because 
'tQ(t)  and  W^(t)  are  known  from  the  initial  conditions  (1.6)  and  (1.7). 


1  zr' 


7IT..1 


Introducing  the  linear  oporator 

0.9)  »-•♦* 

wo  oan  write  the  solution  cf  the  differential  equation  in  a  more 
convenient  form. 


Tho  rooursion  formula  is  then 

(1.6') 


Mv(‘) 


The  function  (  scries)  solving  the  equation  oan  now  be  written 
ao  the  sum  of  the  terms  xvW^,  with  v  oven  (  v  *  o,2,4>...)  and  of 
the  torms  with  v  odd  (v«»  )• 

For  the  ?  with  even  index  v*e  write 

v  +  2  =  2(i 


This  now  summation  index  |i  is  inserted  in  tho  roourBion  for¬ 
mal..  (1.8') 

■  zulin-u  ■  jhyr  B'Vt) m  t^hTs 

Likewise ,  for  odd  v  we  have: 

v+2*2p  +  1 


v+2 


(t) 


■  '  T)U 

2u (2  4+1 )  V*2  4-I 


1 

(24+1)1 


n\(t) 


1 

I2H7T)! 


Hence  follows  the  formal  solution  of  the  telegraphic  equation  (1 . 

°°  *  2V  2^+1  y 

(i.io  »(*.«)  -21  Ttj!DVf(t)  +H.  fsv+iTi D  h(t) 

v*o  va0 

To  obtain  the  oomplete  solution  one  has  to  know  the  corres¬ 
ponding  boundary  and  initial  conditions. 


The  current  j(x,t)  is  obtained  by  integration  of  the  differential 
a cue  cion  (1.1) 

(1.11)  J(x,t)  =  JQ(t)  -  J  (q  *  c  |^)U(x,t)  dx 

'he  proof  that  ( 1 . 1 o )  converges  can  be  given  by  a  majorant  method f 
it  can  be  found  in  / 23/*  p.  H4  et  sequ.). 
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VII. 2  Switch-on  Transients  with  Shorted  Wires 


1.  Initial  and  boundary  conditions 


To  bo  able  to  solve  tho  telegraphic  equation,  i.  e.,to  describe 
the  boundary  and  initial  value  problem  completely,  we  need  the  cor¬ 
responding  boundary  and  initial  conditions. 

We  shall  assume  the  following  conditions: 

Until  time  t=o  there  is  no  current  nor  voltage  in  tho  wires. 

An  alternating  voltage-  U(o,t)  =»  A  cos  wt  +  B  sin  ut  is  applied  at 
^his  moment  t»o.  This  immediately  gives  tho  initial  conditions  for 
the  voltage  funotion  U(x,t)  at  any  x  except  at  x=o  and  for  tho 
current  function  j(x,t). 

(2.1)  U(x,o)  =  o  x  >  o  „ 

(2.2)  J(x,o)«o 

The  initial  condition 

^(x.o)  =9  x  >  o 

follows  from  Eq.  (1.1 ). 

What  we  still  need  are  tho  boundary  conditions  for  both  ends 
of  the  lino  of  length  a.  One  of  them,  for  x=o  (after  an  ordinary 
alternating  voltage  has  been  applied),  is 

(2.3)  U(o , t )  =  A  cos  wt  +  B  sin  wt  »  f(t) 

As  wo  assume  the  line  to  be  sherted,  the  other  boundary  con¬ 
dition  for  x=a  at  the  end  of  the  line  is 

(2.4)  U(a,t)  =  o 


.V 
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2.  Transformation  of  a  fen  expre salons 


Before. introducing  the  initial  and  boundary  conditions  into  the 
formal  solution  for  U(x,t)  (cf.(l.lo))  we  shall  bring  a  few  expres¬ 
sions  into  a  more  convenient  form. 

Wo  apply  the  generalized  Lie-operator  D  of  Eq . ( 1  • 9 )  to  the 
function  f(t)  which  gives  the  variation  of  the  voltage  U(0,t)  at 
the  point  x»o 

D°f(t)  =  f(t)  -  Aqoos  wt  +  BQsin  wt 

D^t)  =  Df(t)  *  A^os  <ot  +  B.jSin  o>t  «= 

2  2 
■  Qa-yw  )A  +  PwBjJ008  wt  +  £-  (3u>Ao  +  (o-yaj  )Bo^sin  yt 

This  gives  the  following  relations  for  the  coefficients  A1  and  B1 
A1  -  (a  -  y  u  )Aq  +  Po)3q 
(a  -  yu?)B  -  pwA 


By  applying  the  operator  D  v  -Times  to  f(t), 


(2.5) 


DVf(t) 


A  cos  u  t  +  B  sin  u>  t 

V  V 


we  obtain  recursion  formulas  for  the  coefficients  Av  and  Bv 
(proof  by  induction)  : 


(2.6)  Ay«  (a  -  yw2)^.-!  + 

(2.7)  Bv  -  (a  -  yo>2)Bv_1  - 


for  v  >  1 


With  the  matrix 


(2.8) 


11  - 


c  -  yw  I 


and  the  corresponding  transposed  (F  (for  the  rules  af  matrix 
calculus  cf.  /24/  )  the  above  formulas  can  be  written  in  matrix 


form: 


-  (1 


v_i 


v-l 


(fiT  )? 


A 


v-2 


v-2 


/A  \ 

All 

l0/ 


"at— 
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After  having  transposed  this  matrix  equation, 

(VBv)  -  (VVL<oVJT  '<Vbo>  flU 

we  insort  the  matrix  (AV»B  )  in  the  first  term  of  the  formal  solutior 
for  U(x,t)  (of.(l.lo)): 


— y,  (A,oos  ut  +  ^  sin  ut) 


V  “  0 


2V 

X 

(2v)! 

Dvf(t)  -  y~f?7) 

v  »0 

2v 

(2.9) 


V  =0 


xf—  X*  ^flts  ,,  .  £1  nv|  (cos  „t\ 

^ - (2v ) !  D  ^  '  v  Ao ’ Bo ' (  2.  (2v)  !  J  \ sin  ut  / 

v*=o  V  v«o  '  '  J 


3.  Statement  for  h(t) 


For  the  function  h(t)  =  irx(o,t)  (see  (1.7))  wo  write 
(2.1c)  h(t)  =  h.(t)  +  h0(t)  »  C  s  ut  +  D  sin  ut  +  hn(t) 

Thoiiwe  substitute  the  expression  (2.1o)  in  (l.lo)  which  gives 

eu  2v  a,  2v+1 

n(x,t)  ,  Dvf(t)  +  ffh,  ♦ 

V  -q  m  © 


(2.11) 


,2v+1 


fern,  tfM*> 


v=o 


In  this  representation, each  term  on  the  right-hand  side  in¬ 
dividually  satisfies  the  telegraphic  equation.  In  the  same  way  as 
•,va  did  with  the  first  term  of  the  above  solution  function  (cf.(2.9)) 
we  can  also  transform  the  second  term: 


008  ut 
sin  ut 


,2v*l  v  /  2V*1  nv  \  / 

<2-12>  L_  T537T)!B  “i(t)  ’  (0o •»o>(X!lfjrTT)iB  )( 

The  first  two  terms  of  Eq.(2.1l)  will  bo  denoted  by  U^(x,t), 
the  last  one  by  v(x,t) 

Tj(x,t)  »  U1  (x,t)  +  v ( x, t )  with 

«  2  v+1 

v(*,t) .  y~  pVy*> 

v=  0 


(2.13) 

(2.14) 


Ho 
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Now,  we  determine  the  coefficients  C  and  D  so  that  the  funotio\ 

o  o 

U^(x,t)  also  satisfies  the  two  boundary  conditions  (2.3)  and  (2.4)* 

One  boundary  condition,  U^(o,t^.  follows  from  (2.11)  when  we 
put  x»o : 

(2.15)  U.j(o,t)«  f(t)  ■  Aocoswt  +BQsinwt 

•  (of. (2. 3)) 

The  other  boundary  condition  is  . 

(2.16)  U^(a,t)  -  0 

The  function  v(x,t)  is  then  zero  at  the  ends  of  the  line  (x=o 
and  x*a), 

Henoe,  we  have  another  boundary  oondition  for  v(x,t),  namely 

(2.17)  v(a,t)  -  0 

Moreovor,  we  must  determine  v(x,t)  so  that  it  satisfies  the 
initial  conditions  (2.1)  and  (2.2)  for  TJ^x.t)  does  not  satisfy  them. 

As  we  shall  see  later,  the  function  U.j(x,t)  constitutes  the 
steady  part  of  the  solution  whereas  v(x,t)  is  the  non-persistent 
part  of  the  voltage  function. 


I 

4.  Calculation  of  the  coefficients  C  and  D 

0  0 


To  make  the  calculation  of  C  onl  D  more  transparent,  we  bring 

c  0 

the  matrices 

00  2v  on  A)  +1  v 

pi  '  ZZxuji  n  *2  “  211  (a, +1)1  0 

V  ■  0  0 

to  their  normal  form  (of.  / 24/,  p.  1 95 ) - 


The  eigenvalues  ^0  matrix 


(2.18)  q1  2  -  (o  -  yj)  t  i.Pu 


are 
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The  matrix  P1  »  cosh  xf?T  oan  he  brought  to  diagonal  form,  be¬ 
cause  the  eigenvalues  q1  and  q arc  different  (0  /  o) 

(2.19)  T"1  cosh  x*y?T  T  -  ■drj7(T“1ftV  T)  -  ^  ■{|jp|(T''10  T)' 


v-o 


**“  .V 


2v 


TSTTi" 

V  -o 

^cosh  xfc 
o 


v-o 


cosh  x\j$2 

After  a  simple  calculation,  we  find  the  matrix 

(2.2°)  T. 


(2.21)  T 


-1  1 


1  1 

-i  i 

!\  i 


-1  ^ad 

and  it 3  inverse  T  -  -j^p 


1  -i 


TtTi  (o’ 


Multiplying  the  matrix  function  (2.19)  on  the  left  by  T  and 
on  the  right  by  T””*  wo  find 


.  '/cosh  +  cosh  xifcT,  i(cosh  xT/q] 

oosh  xF  -  ~  f  T  1  ! 

*  l-i(cosh  xy^j  +  oosh  x|q^),  cosh  x|o]j 


}  +  oosh 

O’  +  cosh  x^J 


After  a  few  transformations 


-  p  -  iq 

yo^=  ^  (a  -  yu  )  +  P«i  -  p  +  iq 


P  iSL  jjtil  +  ljj(a  .  yw2)2  -  (ft,)2  ' 


2, 


*  f^(«  -  wV  -  (ft.)2 


we  find : 

(2.22) 


P1  *  cosh  xjfjf 


Pi 

-qi 


p1, 


oosh  xp.cos  xq,  sinh  xp.sin  xq 
-  sinh  xp.sin  xq,  cosh  xp.cos  xq 


! 


142 
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In  the  same  way  we  also  troat  the  matrix 

1 


(2.23) 


P2  ‘ 


2  v+1  v 

X  .  ov 

T2V7TT! 


v  -o 

The  oulculation  gives  then 


n  2  #  sinh  xfn 


(2.24) 


/  P2  M 

\  -"2  J>2  / 


with 


2^2 

p  +  q 


(  p.sinh  xp.cos  xq  +  q.oosh  xp.sin  xq  ) 


i  ■,  (  p.cosh  xp.sin  xq  -  q.sinh  xp.oos  xq  ) 


2^2 
p  +  q 


With  x*ia,  the  matrix  P1  beoomes 


,2v 


M 


V-0 


(2vJT  °  “  coah  aYn 


/n11  m12  \ 

Vm21  m22  / 


with 


m11  “  m22  “  00sil  aP,oos  act 
m12  "  “  ra2l  “  8inh  aP,ain  acl 


whereas  the  matrix  P.  becomes 

m  _  2 V+1  y 


and 


n 


N  =  _  (I'vtT) !  °  0  2  s  irh  G 


v-o 


11 

l2l 


n 


n 


12 

22, 


with 


n<  „  =  n~„  ■ 


11  ‘22  2,2 
P  +q 


n12  **n21 


2,2 

p  +  q 


(p.sinh  ep.oos  aq  +  q.oosh  ap.sin  aq) 


(p.coai  ap.sin  aq  -  q.sinh  ap.oos  aq) 


C/naidering  the  boundary  condition  (2.1*)  we  obtain 

Ui  (a,t )  -  (A  ,B  )  M /C?SW!\+  (C  ,D  )  n/°°S  - 

1 v  ’  '  v  o*  o'  l  sin  »t  f  '  o  o'  |  sin  wt  J 

-  (Vll  +  Bo*21  +  °on11  +  Don21)oos  *  + 

*  <Vl2  +  BoD22  *  °on12  +  Bon22)"ln  "*  ‘  0 

The  coefficients  of  cos  wt  and  sinw  t  must  vanish  in  order  that 
the  above  relation  is  satisfied  at  any  instant  of  time  t.  The 
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coefficients  Cq  and  DQ  arc  then  readily  found  by  Cramer's  rule. 


.1 . 2  •?  -  C  - 
o 


Ap(m11nr.  m1 2n2 1  ^  +  ®0^Vl"22“  r^22n2 1  ^ 


n12n21  “  n11n22 


(r.,25)  D0  « 


A0(n1 ln12*  ra12n11>  +  Bo^m21n12“  m22n11 > 
nlln22  "  n12n2l 


7>.  Calculation  of  the  voltage  part  v(x,t) 


The  oend.it  ion 


cj  _  2y+1 

\ -  a  _y 


,^•27)  v(a.t)  =  XZ  72—).  h2(t)  "  0  (of. (2.17)) 

cc.n  ho  fulfilled  when  the  function  of  time,  DVh2(t),  is  formally 

nr.  turned  as 


2  2  ~v 

(2.28)  DVh2k(t)  -  (..1)V.  h2k(t)  .  pS-S-j 


for  k  *  1j2y5>*** 


Zhir  relation  must  bo  inserted  in  v(x,t): 

*£l  7S7T7,  ®V<‘> 


V  r-.O 


v-o 


Using  the  series  expansion  of  the  sine  and  superposing  v(x,t) 
v;o  find 


(2.29)  v  ( x  j  t ) 


<*.*)  -K  fe  W*) 


k-*1 


knx 

a 


The  Relation  (2.28)  is  equivalent  to  the  equation 


(2.;o)  Dh2k(t)  -  -  h9Jt). 


.2  2 
k  n 


2k’ 


’A hone c  one  can  determine  the  functions  h2k< 


.2  2 
k  n 


-■aeJX)  "  c:h2k(t)  +  P  *5t  “jkO  +  *  -T  h2fc(*>  -  -  h2k<*>  2 

0t  a 


V 
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.2  2 

yhg'k^)  +  Ph2k^  +  +  — T~)  h2k^^  “  0 

a 

This  homogeneous  differential  equation  has  the  solution 
-  ££ 

(2.31)  h2k(t)  -  e  2y  (Ckcos  wkt  +  D^sin  w^t) 

\  •V(“  +  H')f -"J- 

(the  case  w2<  0  oan  easily  be  inolused  in  the  oaloulation). 

The  funotion  hg(t)  can  also  be  obtained  by  superposing  the  h2t(t; 


(5.32)  h2(t) 

Inserting  (2.32)  into  (2.29)  we  obtain 

„ii  . 

(2.33)  v(xft)  -  e  2y  2^1  H  (Ck°0B  V  + 


Dksin  Wj^t )  sin  — 


6.  The  solution  U(x,t) 


Inserting  (2.33)  into  (2.13)  we  obtain  the  solution  function 
(2.34)  U(x,t) 


<w(!;  X)  fe:i) +  ^{X  4 


+  e 


.1 1 

2y 


k-1 


r—  (  C.  cos  w. 
kn  \  k  l 


t  +  D^sin  t 


\  i  knx 
) sin  — 


The  coefficients  C,  and  D,  follow  from  the  initial  conditions  (2.1) 

k  k 

and  (2.2). 

We  insert  the  condition  (2.1)  in  (2.34)  and  we  obtain: 


U(x,o) 


2  42W1 

q2  P2 


k-1 


a  _  .  kwx 

—  C,  sin  — — -  »  o 
kn  k  a 
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The  C.  can  be  regarded  as  the  Fourier  coefficients  of  the 

Ki*  K  f  tffY  I  w  is  1# y 

orthogonal  system  sin  »  ^cr  >  C^sin  — —  represents  a 

•n _ J _ _ t _ _ XI _ i1.. k35  1 


Fourier  series  for  the  function 


{■(w(p.;  :;)f 


>.  -  (o-’-> 


Hence, we  obtain  the  Fourier  coeffieients 


<2-”>  57  Ck  “  ‘  a  J{(io’Bo>  ( -ij  +  (0o’Do)  (-0  }  8i"  ^  ** 

o'-  v  *  '  ^ 


The  result  for  the  coefficients  Ck  is 


(2.36)  C. 


J  ui 


(x,o)  sin  ,  <jx 
a 


To  determine  the  coefficients  we  have  to  consider  the  initial 
oondition  (2.2). 

U+(x,o)  =  Ui+(x,o)  +  v  (x,o)  *  o 


Calculation  gives 

V*'°>  -  (vvf-J,  pj)(  i)“  +  (c 


o^(X  P22)(?)“- 


k»  1  ’  • 

,  2  , „  _  n  a  /o\  knx 

rz  <w  —  ( i)  ukain—  ■ 

I 


Here,  2L_  \  *7  ^ 
series  of  the  function 


rj  J/«r  Y 

I>k  —  sin  — can  a<>ain  be  assumed  to  be  the  Fourier 


{-  (Ao-v(p ’/•“  -  <*.•». >(,*)• ♦ 

+  f?  [-(V-.)  (!;,)  -  <co*Do>  (!q?)  _ 


VJith  the  Fourier  coefficients 


<2^>  57  Vk 


f  j'-|  -  -  |y  'J,  (x,ajBin 


1 40 
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7.  The  solution  J(x,t) 


Integrating  (1.1)  we  obtain  the  solution  for  the  ourrent  j(x,t) 
(2.39)  J(x,t)  -  -  J  (q  +  o  -|^)  U(x,t)  d x  +  JQ(t) 

With  the  abbreviations 

pi(x)'(-q1  pj)'ooshxVn 

<\  _L 
A  -ft  2sinh  xVn 

Pj(x)  »  Jp2(x)  dx  -  ft  cosh  xW 
vo  find;  .  .  f  . 

<*•*•>  J(M)  ■  -  S{(A0>B0)  lazily  (C0.»e)  P5p  3)- 

-I  k„* 1 

033  —  f  - 

— 1  J 


P2U) 


(■ 


<>2  P2 


_  M 

2  y  /a_ 

k-1 


C,.cos  c^t  +  D^sinuvt  |cos 


-  °{<W  F2  “  +  «vV  P5  ('co"  s)* 

-  a  .  *  \  ?r  — 

\  tJc' 


+  f-  .'  2» 

2y 


C,  cos  (u^t  +  D^sini 


kxx 

a 


-  e 


-|i 

2y 


!(tn]f  [-  Ck‘“k3in  V  +  Vk  003  v]003  +  Jo(t) 


The  integration  constant  J  (t)  can  be  found  from  the  relation 

o 


for  j(o,t)  which  can  be  derived  direotly  from  the  basio  equation  (l.‘. 

The  time  function  j(o,t)  satisfies  the  differential  equation  ('.') 
(1.1)  U  (o,t)  -  h(t)  -  -  (r  +  1  l^-)  J(o,t) 


lJt(o,t)  +  rj(oft)  a  h(t) 


(2.41) 
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Its  solution  is 


C  1 


(?.42)  J(o,t)  »  -  -2 - —  (  £.  coa  ut  +  u  sin  ut)  - 

r  +o)l  1 


H  1 
o 


-  1  g  g  (  T  s  in  u)  t  -  o)  c  o  s  o)  t )  - 

T  +  wl 


C 


-5^7 


.  £1 


1: 


k-1 


<?-&>*  -4 

-  it 

2y 

e  _ 

(  £  .  JL)2  +  o)2 

K  1  2y]  k 


"  L('i' "  $r-  008  “k*  +  Vin“kl  J  - 


I  (f  -  ly)  ‘  V00  “/I  +  Ke 


The  constant  K  results  from  iho  initial  condition  j(x,o)  =  0 

(of. (2.2)) 

(2.43)  J(o,o)  -  o 

C.r-D.lu)  +  Ck  (y  -  |y)  -  Dku»k 

k=1  - 17)2  +  “k  ] 


(2.44) 


o  o 
“  "  2 


r  + 


2i8 


Hence,  ne  find  the  original  integration  constant  J  (t) 


(2.45)  J„(t)  -  J(f  ,0  *  »|(C0.B0)  a_YaSnSt)  ' 


fit 


-  e 


+  0 


)ksin  ^ 


a} 


0t 

+  JLe" 2? 

2  Y 


-  M 


-  e 


Fr^llv08  «i# + vin“k‘3  - 

_  (  -  C^w  ^sin  u jjCos 


V 


8.  Numerical  examples 


Two  numerical  examples  shall  illustrate  the  switoh-on  transient 
problem  in  a  shortecl^We  have  examined  the  maximum  values  of  the 
power  function  at  the  point  x*o  in  dependence  on  the  phaso  angle  t 
of  the  applied  voltage 

U(o.t)  •  u  cosf  « t  -  t)  ■  A  cos  tii). t  +  B  sin  ^t 
'  f  '  0  o 

during  the  transient  process.  Calculations  ^ere  performed  at  the 
ZUSE  Z23  computer  of  Innsbruok  University. 

In  the  first  numerical  example ,  the  electrical  constants  per 
unit  lenght  were  chosen  as  follows : 

resistance  r  «  1 

induktive  reactance  1  -  o,2 

capaoitive  reac¬ 
tance  o  -  o,oo2 

leakage  q  «  o 

This  result  (  see  table  )  3how3  that  there  is  a  notable  strong 
resonance  and  superposition  effekt.  For  t  -  12o°,  the  power 
funotion  main  peaks  at  x^o  increases  to  more  than  2o5$  of  the 
maximum  value  in  the  steady  final  state. 


lenght  of  wires  a  -  1 

angular  frequency 
of  voltage  u>  -  3oo 


+) 


The  corresponding  programs  of  the  numerical  examplos  are 
contained  in(/45/p»  59  ff.) 
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The  following  toble  ahowa  the  result 


T 

A  eloo.oos  T 

0 

B  sloo.sin 
0 

0 

1  00 

0 

l0 

98,4800772 

17,3648176 

2o 

93,9692619 

34,2o2o142 

3o 

86 ,6 o254o 1 

50,0000000 

40 

76,6o44441 

64,2787608 

5o 

64,2787608 

76,6o44441 

60 

50,0000000 

86 , 6o2 540I 

j  7o 

34,2o2o142 

93,9692619 

80 

17,3648176 

9a, 4800772 

9o 

0 

loo 

1 00 

-17,3648176 

9fi,48oo772 

11o 

-34,2o2o142 

93,9692619 

12o 

-5o,ooooooo 

86,6o254oi 

1  50 

~64,27876o8 

76,6o44441 

1 40 

-76,6o44441 

64,27876o8 

15o 

-86,6o?54o1 

5o , 0000000 

1  60 

-93,9692619 

34,2o2o-'42 

17o 

-98,48oo772 

17,3648176 

I80 

-loo 

i 

0 

19o 

2oo 

2lo 


90,48oo772 

93,9692619 

86,6o254o1 


-17,3648176 
-34»2o2o1 42 
-5o ,0000000 


Main  peaks  of 
power  function 
during  transient 


Peaks  in 

final 

state 
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Fig.  Peaks  of  the  ,-oser  function  N(o,t)  ■  j(o,t).u(o,t)  versus 
phase  angle  t  of  the  applied  alternating  voltage. 


Max  N(o,t) 


4ooo 

3ooo 

2ooo 

looo 


transient 


steady  final  state 


o  3o  60  9o  12o  15o  18o 


phase  t 


The  following  values  sere  chosen  in  the  second  numerioal  example : 
resistance  r  ■  1  lenght  of  sires  a  ■  1 


resistance 
inductive  reactanoe 
capacitive  reactanoe 
leakage 


1 

1 

o,o1 

o,o1 


angular  frequency 
of  voltage 


With  N(o,t)  tabulated,  the  following  result  sas  obtained  for  the 
peaks  of  the  poser  function  N(x, t)  at  the  point  Xr-o 


Phase  t 


I 

V  ■loo.cos  T 

0 

B  =1oo.sin  t 

0 

Main  peaks 

during 

transient 

Peaks 

final 

loo 

0 

1144 

928 

98,48oo772 

17,364817 6 

1166 

928 

95,9692619 

54,2o2oH2 

1127 

928 

86 ,6o254o1 

50,0000000 

1o28 

928 

76,6o44441 

64,2787600 

1o29 

928 

64,27876o8 

76,6044441 

1o37 

928 

50,0000000 

86,6o254o1 

1o4o 

928 

34,2o2o1 42 

93,9692619 

1o42 

928 

17,3640176 

98,43oo772 

I068 

928 
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Phase  t 

A  -lOOeCOS  T 

0 

33  “loo.  irin  t 

0 

Main  peaks 

during 

transient 

Peaks  in  steai; ■ 
final  state 

1  8o 

17,3648176 

98, 4800772 

I068 

928 

9o 

c 

1  oc 

1o43 

928 

1 

1  OQ 

-17.3648176 

98,48oo7?2 

1o42 

928  | 

1  1  0 

-34 ,2o2oH2 

93,9692619 

lo38 

920 

12o 

-50,0000000 

86,6o254ol 

I036 

928  | 

13o 

-6/.  ,2737608 

7S,6o4444l 

lo12 

928  | 

.  Ho 

-76,6044441 

64,2787608 

1o24 

920 

! 

-86,6o2 54o1 

50,0000000 

1o27 

928 

j  uO 

-93,9692619 

34,2o2o142 

I008 

928 

1  |7o 

.•98, 4300772 

17,3648176 

1o36 

928 

1  i2o 

1 

-Ico 

0 

1144 

926 

|  13o 

-98, 48oo772 

-17,3648176 

1166 

928 

:  2oc 

-93,9692619 

-34 ,2o2o1 42 

1127 

928 

rhi&  example  shows  a  clear  dependence  of  the  main  peaks  on 
i)Lftpp  angle  of  the  applied  alternating  voltage. 


Tig.  Main  peaks  of  the  power  function  H(o,t)  during  the  transient 


an  1  peaks  of  N(o,t)  in  the  steady  final  state  versus  phas 
a>agle  t  of  the  applied  alternating  voltage 
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VII. 3. 


VII. 3  Switoh-on  Transient 8  rith  Open  Wires 


The  transient  processes  with  the  wire  ends  open  oan  be  treated 
mathematically  analogous  to  the  case  of  shorted  ends. 


1.  Initial  and  boundary  conditions 

We  shall  assume  the  initial  oonditione 

(3.1)  U(x,o)  «=  o  x>  o 

(3.2)  J(x,o)  -  o 

(5.3)  Jt(*,o)  -  o 

and  the  boundary  conditions 

(3.4)  U(o,t)  -  A  cos  ut  ♦  B  sin  ut 

(3.5)  J(a,t)  -  0 


Liko  in  (1.1  e),  the  formal  oolution  for  the  current  is 

c2v+1 

[57fT)! 


2v  go  2v+1 

TCTf  D^Tt)  +  75^*  ^h(t) 

'  '  v»0 


(3.6) 

00 

where 

(3.7) 

J(x.t)  -  5 

V**C 

f(t)  -  J(o,t) 

and 

(3.8) 

¥(t)  -  Jx(o,t) 

The 

function  'E’(t)  0 

and  from  (3.4) 

(3.9) 

F(t)  -  Jx(o,t) 

=  X  cos  mt  +  U  sin  wt 

whore  wo  have  put 

-Bwo  -  $A  «  S’  and 

Au>c  -  qB  *  IT 
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2.  Statement  for  j(o,t) 


Fox  the  funotion  j(o,t)  -fie  set  formally  (of.  VII.2.3) 

(3.10)  j(o,t)  -  7(t)  =  ^(t)  +  f2(t)  - 

=  C  oos  u  t  +"D  sinwt  +  71  (t) 
002 

which  we  insert  into  (3*6) 

(3.11)  J(x,t)  -  72vT!D^1^t^  +  TSv) + 

v  »  0  V  «  o 

co  2v  +1 

\»=  0  '  • 

We  call  the  sum  of  the  first  and  third  terms  J.,(x,t)  and  tran3 
form  them  in  the  same  way  as  in  (2»?) 


(3.13) 

(3.11‘) 


Hr.rtri  in  (3.11)  wo  write 


2m 


w(x,t)  -  >  ‘  D  f,(t) 

v  ■  0  '  v  '  ' 

j(xft)  a  J^Xft)  +  \v(x,t) 


The  coefficients  TT  and  7  must  be  determined  so  that  the 

c  0 

function  J^(xft)  satisfies  the  boundary  condition  (3*5) 

(3.14)  J1 (a>t)  -  o 

J^Xjt)  does  not  satisfy  the  initial  conditions. 

The  function  w(x,t)  must  bo  determined  so  that  the  initial 
conditions  are  fulfilled  and  tha.  w(x,t)  becomes  zero  at  the  end 
of  the  line  (x=q). 


/ 


1*5/! 


VII. 3. 


Thus,  wo  have  fux'lliHT*  bound  ary  condition  for  w(y,t) 
' 5)  w(n,t)  =  0 


3.  Calculation  of  coefficients  0  ,  D  and  of  tho  function  w(y,t) 

o  o 


Using  the  boundary  conditions  (3.14)  we  obtain  (of.  VII. 2. 4) 

i 

t 


(ms)  »,(.,«).  :;*)  (s : 


:) 


003  ut 

sin  ut 


).° 


This  relation  holds  at  any  time  t.  Therefore,  the  coefficients 

of  sin  ut  and  cos  u  t  must  be  zero.  The  coefficients  A  and  B 

arc  given  by  the  boundary  condition  (3 • 4 ) *  whence  the  coefficients 

(T  and  IT  can  bo  found  by  means  of  Cramer's  rule:  the  result  is 
oo  * 

,  ,  r  r.(|l1lV  -  n12°Vl>  *  -  n22m21> 

o  ra12m21  ”m1 1m22 

— .  Vn1 lm1 2  ”n12ra11^  +  ^o^n21m12  "  n22n1 1 ^ 

Do - 


22  11  21  12 


The  condition  (3*15) 

w(a,t)  =211  72UT»  a  0 

V  S3  0  '  ' 

is  satisfied  if  we  put 

(3.18)  Dvra(t) .  (-i)v[(2fe±1)->]a’r21[ 


for  k  *  o  ,1 ,2 , . . . 


V') a  tl  f^T'  (  if')  •  f2k(t) 


v=  o 


cos 


(**&  «)  nk(‘) 
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Superposing  too  -j,  vx,l;)  the  current  part  \v(x,t) 


(5.19)  v/(x,t) 


k=o 


(2k+l) 


2a 


itx 


Tho  relation  (5 .13)  can  bo  written  no  a  differential  equation 
of  the  form 

(3.20)  D^j^t)  -  3f^(t)  P  •,aV*’f9lr(t)  +  I  f ?!,(*)  = 


2k' 


dt  2k' 


at 


■  (-1>-~2k(^2T'  *) 

^2"  (t)  +  pr2k(t)  *  [ "  ]  r2k(t>  ■ 0 

The  funotions  f2k(t)  can  then  be  found  from  this  equation. 


(5.*1)  f21;(t)  -  e 

-.iih 


» £i 

2  V 


(S^Oos  w^t  +U,c3in  w^t) 
/2k+1  \2 1  1  g2  1 

I  +  ^--Vj7  -^r 


,h.j- 


v  ha’c  nnuumed  that 
+  )>$ 


A>  The  Point  ion  j(x»t)  and  •!  (:< ,  t ) 


Inserting  (3.21)  in  (3.1?)  wo  obtain 
(3.22)  j(x,t)  =  J1(x,t)  +  W  (  X  »  t  )  = 


2i/coa  a 


-  JS. 


2^  ^  °  (7,  COG  u>kt  +  Uj^sin  U^tjcos  ^  rj; 


k=o 


+  0 


I 
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Thj  coefficients  and  <’an  be  onlculatcd  by  the  sane  ronaonin,' 
an  in  VII. 2. 6. 

Considering  the  initial  conditions  (3«2)  and  (3*3)  we  obtain 
the  result 

(3.23)  7k  -  -  f  J  J^Xjo)  cos  ■  dx 

0 

\  ■  -  ^  j{Ji  +  V, (*.»)}«>•  t**#* 

0 

The  solution  for  the  voltage  U(x,t)  is  found  by  integration 
of  (1,2): 

(3.25)  U(x,t)  -  -  j"  (rj  +  1  |^)ax  +  0o(t)  . 

.  - rf(c,5)  p ,ha»t)+  a  j) p,  hs )  + 

i  '  o’  o'  2\sin  yt  j  o’  o'  3  \sin  uit  / 

■^Tlx  <V08  “k*  +  \8in  "k*)  8in  - 


+  o 


6t 

••  — ■  ••  no 

2y  v~ 


k*o 


-  l{<W  P2  ('«”  lit)  “  +  <VBo>  P3  (‘cosl  J  )  “  * 

-  2j  °  2r  ZI  (V08  +  Tkaln  V>  (2ktl)«  8ln 

-  2~-  X 

1  2T.'  rji^fyi  (  -  ^"k  8ln  “k*  +\“k  ooo  ntO. 

k-o 


+  0 


.sin 


(2k+1  )nx 


2a 


}  +  Po 


(t) 


•Jherc 


Uo(t)  -  U(o,t)  +  r  {x}  ♦>{*} 


A  cos  ut  +  B  sin  ut  + 
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